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INTRODUCTION 


The flow field within advanced axial flow turbomachinery for aircraft pro- 
pulsion systems, specifically high bypass ratio turbofan engines, is character- 
ized by the presence of mixed subsonic, transonic and/or supersonic regions, 
multiple shock waves, shock-boundary layer interactions, and significant ef- 
fects of three-dimensionality and unsteadiness of the flow. These aspects of 
the flow field are not easily accommodated in the conventional analytical and 
numerical methods available for application to turbomachinery, and the design 
and analysis of advanced systems has been correspondingly constrained. A fairly 
recent review of methods for steady cascade analysis (i.e., two-dimensional 
compressible flow) is presented in Reference (l). In the area of transonic and 
mixed subsonic-supersonic flows, the most promising approaches appear to be the 
relaxation method for solution of the transonic small perturbation equation de- 
veloped in Reference (2) and finite-difference solution of the complete system 
of equations for an inviscid compressible flow, as applied to a cascade in 
Reference (3), for example. A steady three-dimensional transonic solution for 
a blade row obtained by the finite-difference method is presented in Reference 
(k) , and a corresponding solution for a wing using the relaxation method is 
presented in Reference (5). Analyses of unsteady flows in turbomach i nery have 
been limited to solutions of the small disturbance equations, usually with a 
view toward description of the acoustic output of turbomachinery; cf. References 
(6), (7) and (8). Extens ion of the relaxation method to solution of the tran- 
sonic small disturbance equation for an oscillating airfoil is presented in Ref- 
erence (9)- Thus, while important advances have been made in the last few 
years, a considerable amount of work remains to attain a comprehensive predic- 
tive capability in the subject problem area. 

% 

The present study has been directed toward development of a numerical method 
of solution of the complete unsteady equations of motion for a compressible, two- 
dimensional flow through a blade row (either rotor or stator) or a stage (both 
rotor and stator) of a compressor or fan. The objective is attainment of a 
steady solution for a single rotor or stator blade row or a periodic solution 
for an interacting pair of blade rows In a stage. Either case may include 
mixed subsonic, transonic and/or supersonic flow containing embedded shock waves. 
The initial effort is described in Reference ( 1 0) . 
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The analysis is formulated with respect to a blade-to-blade stream sur- 
face, as depicted in Figure (1) for one passage of a blade row. The inlet 
and discharge boundary conditions are applied at axial stations some distance 
upstream and downstream of the blade rows. The selected boundary conditions 
assume subsonic axial velocity at both stations, but admit either choked or 
unchoked operation of the blade row or stage. In the case of a single blade 
row, or the equivalent infinite cascade on a b lade-to-blade stream surface, 
it is immediately evident that the steady solution must possess a periodicity 
in the circumferential direction with a fundamental period of 2 tt/N^, where 
is the number of blades in the row. (it is implicit that the inlet and dis- 
charge boundary conditions also admit this periodicity condition.) If the 
blade row is rotating these considerations also pertain in the rotating frame 
of reference. Thus the computat ional domain need only encompass that fraction 
of the flow annulus containing the flow through a single blade-to-blade passage. 
The locations of the boundaries of the blade-to-blade passages may be defined 
arbitrarily so long as their spacing corresponds to the blade pitch. Down- 
stream of the blade row the blade slipstreams represent natural boundaries of 
the b lade-to-blade passages, since their spacing is identically the blade 
pitch, and certain components of the solution may be discontinuous across the 
slipstreams. Upstream of the blade row the boundaries have been conveniently 
defined as projections of the mean camber lines. 

In the case of a stage consisting of a rotating blade row and a stationary 
blade row, a set of blade-to-blade passages may be defined for each blade row 
In accord with the above considerations, but the noted circumferential peri- 
odicity condition will not apply. As has been established from considerations 
of the acoustic problem, cf. Reference (11), the flow pattern will rotate with 

an angular velocity of N D fi/(N -N_) and have a circumferential period of 

K SR 

2tt/(N s -N r ) . Only in the limiting case of N R =N S> which is usually avoided in 
practice, will identical solutions occur in each b lade-to-b lade passage at 
any instant. In the typical stage (N g >N R ) the flow pattern will rotate in the 
opposite direction of the rotor. The solution in any particular passage at 
one instant can, therefore, be related to that in another passage at an earlier 
time. This phase lag forms the basis of a cyclic procedure developed herein 
for relating the conditions along the boundaries of the computational domain 
(which is composed of a set of blade-to-blade passages for the two blade rows) 
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to the solution within the domain at an earlier time. The cyclic procedure 
will be described in detail in a following section. 

The total pressure losses due to boundary layers on the blades will be 
carried in the wakes of the blades. Thus, a significant contribution to the 
unsteady aerodynamic interaction between blade rows and resulting acoustic 
signals may be attributable to passage of one row through the viscous wakes 
of the other row. Accordingly, an approximate representation of the boundary 
layers on the blade surfaces and the blade wakes is also incorporated in the 
present model. The considered viscous effects include reduction of the blade- 
to-blade passage area due to the boundary layer displacement thickness and 
turbulent diffusion of the corresponding momentum defect (i.e., total pressure 
loss) in the wake. 


INVISCID FLOW ANALYSIS 
Fundamental System of Equations 


The basic system of equations from which the present flow representation 
derives consists of the statements of conservation of mass; 


9p 

at 


+ 


V-pV 


0 


conservation of momentum; 


[W Vp 

Dt p 


0 


( 1 ) 


(2) 


and conservation of energy; 


De 

Dt 


D ,U 

p w ( j } = 


where : 


_D_ 

Dt 


_9 

at 


= e- + v-v 


Assumption of an inviscid gas is Implied by the absence of shear stress and 
heat conduction terms in Equations (2) and (3). In addition, a thermally 
and calorically perfect gas is assumed, for which: 


(3) 

(V 


h 



p - pRT (5) 


e “ IFFTTF 


(6) 

The energy equation can be stated in a somewhat more convenient form through 
combination of Equations (l), (2) and (3); 


+ V-pfa 
at 

= 0 

(7) 

where: 

E = e + 1 

v 2 

(8) 

H " h + I 

V 2 E 

(9) 

h = e + p/p 

do) 

Furthermore, Equation (3) 

can also be written as: 


fif. = 0 

Dt 


(11) 


where: 

p = kp Y exp(S/C y ) (12) 


An important distinction between Equations (7) and (3) or (11) should, however, 
be recognized. Equation (7) is in divergency or "conservation law" form, 
whereas Equations (3) and (11) are in a non-conservative form. The divergency 
or "conservation law" form of the equation may be applied across a shock wave, 
but the non-conservati ve form can not. 

The only essential simplifications to this system of equations which are 
introduced in the following development of the flow model are a reduction in 
spatial dimensions of the problem from the general three-dimensional form 
used above to a particular two-dimensional form, and use of small disturbance 
approximations in the acoustic far-field of the machine. Derivation of the two- 
dimensional system is outlined below; discussion of the acoustic solution is 
presented in Volume 111. 
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Absolute Stream Surface Equations 


Consider a flow annulus as sketched in Figure (2). The curvilinear dis- 
tance along the intersection of the mid-line of the annulus with a meridional 
plane is denoted by m, termed the meridional distance. The distance normal to 
the mid-line in a meridional plane is denoted by n. The circumferential coordi- 
nate 0 is considered positive in the coun te r-c I ockwi se direction when viewed 
down the positive z axis. The thickness of the annulus b is assumed to be 
small compared to the radius r; hence the n component of the velocity vector 
and all variations in the n direction are neglected. Accordingly, the annulus 
is termed a stream surface. Transformation of a system of equations of the 
above form to the considered two-dimensional coordinate system is outlined in 
Reference (12). Application of the some procedures to Equations (1), (2) and 
(7) yields: 


•ip 

?)t 


+ 


>i p V 


m 


3m 


+ 


, ^B 

r 36 


+ 


pV , , 
m d rb 

rb dm 


n 


3V 

m 

at 


+ 


v 

m 


3V 

—01 + 

am 



r 


3V , . „2 

m + I iP _ Vh d_r 

30 n 3m r c j m 


av 


6 


at 


+ 


v 

m 



av. 


r ae 


J_ 3£ 

P r 36 


W dr 
r dm 


apE + _3_ 

3 1 3m 


(r V H) 
m 


+ 


1 JL 

r 36 


(pV fi H) 


+ 


pV H 
m 

rb 


d rb 
dm 


0 


(13) 

(14) 


(15) 

(16) 


Equation (ll) transforms to: 


3S 

at 


+ 


v 

m 


as 

3m 


+ 



r 36 


n 


07) 


Equations (13) and (16) are in conservation form. Equations (14) and (15) 
can also be cast in this form by multiplying them by p and substituting Equa- 
t i on (13): 


3pV m 

01 + 

at 


3 

3m 


p» +7&<»W 


pV 2 . . 
^ m d rb 

rb dm 


d r 
dm 


(18) 
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V v , 

ni 0 dr 


09) 


a P v e 

at 


+ (pv v r ) 

dm m 6 


, U h/ 2 
+ 7 7e {oV e 


p) 


pVV 
m 6 

rb 


drb 

dm 


P 


r dm 


Thus, the statements of conservation of mass, momentum and energy given 
by Equations (13)* 08), (19) and 06) constitute the conservative form of 
the governing system of equations, and Equations (13), 0*0, (15) and (17) 
represent the non-conservat i ve form of the same system of equations. The con- 
servative form will be used in the interior of the compu ta t iona 1 domain, 
whereas the non-conservat i ve form will prove to be more convenient to employ 
at the boundaries. 


Relative Stream Surface Equations 


In the analysis of the flow through a rotating blade row it is advantageous 
to express the governing equations in a relative coordinate system which rotates 
with the blades. Therefore, the following additional coordinate transformation 
is i n t roduced : 


t - t (20) 

x = m (21) 

y = r (n-St) at x=constant (22) 

and the velocity components in the x,y,t system are correspondingly defined: 


u = V 

m 

v = - i'ir 

A relative total enthalpy and relative total energy are defined as: 


(23) 

<2M 


H 


H - 


\i rV, 


(25) 


-Note that the "relative total enthalpy" defined herein is sometimes referred 
to as "rothalpy" since it is not the total enthalpy which would be measured 
in the rotating frame of reference. However, it is a quantity which is con- 
served along streamlines in a steady rotating flow, and in this respect it is 
analogous to the conventional total enthalpy in an absolute frame. 
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E 


H 


p/p 


( 26 ) 


The following system of equations, expressed in both conservative and non- 
conservative forms is thereby obtained: 

Continuity 


3p dpu 9pv 

3t 3x 3y 


pu drb 
rb dx 


St Teamwise Momentum 

|£H + 3(pu 2 tP) + 3£i^ = . pui drb + p(v+s2r) 2 idr 
3t 3x 3y rb dx r j x 


or, 


3u 

3t 


+ 




I iP 

p dx 


(v+Qr) 


j_ d_r 
r dx 


Circumferential (Angular) Momentum 


(27) 


(28a) 


(28b) 


2 

3pv , 3puv , 3 (pv +p) 

at - IT 


puv 

rb 


~~ - pu(v+2Qr) _L ijl 
dx r ox 


(23a) 


or. 





1 iE 

P 3y 


u(v+2S2r) 1 di 


(29b) 


Energy 

3pE 3puH 3pvH 

3t 3x 3y 


puH drb 
rb dx 


(30a) 


or, 


as 

at 


+ 




(30b) 


It Is pointed out parenthetically that the above system of equations can be 
applied in the relative (rotating) frame as stated, or in an absolute (sta- 
tionary) frame by setting 0=0 and dropping the prime superscript on H and E. 
In addition, the standard two-dimensional equations of motion are recovered 
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the rotating cascade. It is assumed in the formulation that the second row 
always nas an equal or great*- number o' blades compared to the first row, this 
conforms to standard compressor design procedure. The frame of reference of 
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cludes within it the lower surface of the blade. Similarly, bo refers to the 
1 owe r boundary of the doma.i , wn-ch include:- the upper blade surface. 
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FIGURE 3. BLADE-TO-BLADE COORDINATE SYSTEM AND GRID NETWORK 




The lateral boundaries of domains 4, 6 and 7 are the instantaneous loca- 
tions of the blade slipstreams, 8 (m,t). Each domain is mapped from its shape 
in physical space into a unit square by defining stretched meridional and cir- 
cumferential coordinates, a and v, given (in an absolute reference frame) by: 


a = 


m-m. 

i 


m. , . -m. 
i+l r 


( i n the i ^ doma i n) 


(3D 


9-0 


8 -e„ 

u £ 


where m. refers to the location of the upstream boundary of the i 


th 


of the upstream boundary of domain 1 , , and the location of the downstream 

boundary of domain 7, m^, can be selected arbitrarily. However, 
and are necessarily the planes of the leading and trailing edges of the 
first and second blade rows. The locations and m^ are defined to lie one 
chord length upstream and downstream of the first and second blade rows, re- 
s pec t i ve 1 y : 


(32) 


domain, as 

measured on the stream surface from an arbitrary axial station. The location 


(33) 
(3M 

If and are selected as the inlet and discharge stations at which the 
boundary conditions are to be applied, domains 1 and 7 are not used. Other- 
wise m^ and mg represent the inlet and discharge stations, respectively. 


m 2 = 2m 3 


m = 2m, 

/ b 


Transformation of Variables to Computational Domains 


The relative upper and lower boundaries of the domains, y^ and y^ , are 
functions of meridional distance, x, and may be functions of time, as v/e 1 1 , In 
domains k, 6 and 7 where they represent the Instantaneous slipstream contours. 
Thus, the final transformat ion from physical space (x,y,t) (relative or absolute) 
to computat ional space (a,v,T) is defined by: 


T 


a 


t a /L 
o 


x-x . 

i 



(35) 

(36) 
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(39) 


(40) 


where : 


y = v 

1 T 


9 y u 




(x i+r x i> 


= v 


3y u + n i 3y * 

- + (1"V) 


8ct 


3a 


(41) 

(42) 

(43) 

(44) 


and : 


x. = m. 

i i 


y = r 6 at x=constant 

' u u 


y^ = r at x=constant 


(45) 

(46) 

(47) 


L I s an arbitrary reference iength and is a reference speed of sound. In 

addition, a reference pressure P q will be introduced below to complete the 

non-dimensional izat ion of variables. Particular values for L, a and p will 

' o o 

be assigned later. It should be noted that a and v are non-or thogona 1 coordi- 
nates. The transformation is non-singular for x y ^ 0. 

3 a'v 


The correspondence between the physical and computational domains is indi 
cated schematically in Figure (4). Within the computational domain Equations 
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FIGURE 4. SCHEMATIC REPRESENTATION OF PHYSICAL AND COMPUTATIONAL DOMAINS 




( 27 ), (28a), (29a) and (30a) are expressed concisely by: 


9e 

9t 


<*£ * -U '♦ + *& ♦ h 


where 


e = p x 


1 

u/a 

o 

v/a 

o 

E"/a ; 


ua 

f = o — c 


1 


(u+p/pu) /a 


v/a 
H /a 


va 


g = p 


(u/a o> 

(v+p/p v) / a 


H /a 1 


m 


( 49 ) 


( 50 ) 


( 51 ) 


and 


h 


A 

B 


Pua L 
o 


P 


o 


/I drb\ 
' rb dx } 


x 


1 



0 

u/a 

o 

. a 2 L 
t P dr o 


( v+fir) 2 /a 2 
o 

v/a 

o 

r dx p 

o 

v 

X 

-u(v+2ftr)/a 2 

Q 

H /a 2 
o 



0 


( 52 ) 



( 53 ) 

( 54 ) 
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c = 


(55) 


D - ±- (56) 

Note that f=f(e), g=g(e) and h=h (e , r , b ,ft) thus f, g, and h are known functions 
of the basic dependent variables contained in e. 

The above form of the governing equations is used at all interior grid 
points, as will be outlined in the following section. The non-conservative 
form, i.e., Equations (27), (28b), (29b) and (30b), is employed at the boundary 
points in local orthogonal coordinate systems which will be discussed in con- 
nection with the boundary conditions and boundary point solution algorithms. 

Interior Point Solution Algorithm 

Each of the seven (7) computational domains is spanned by a rectangular 


grid network, having a mesh size of Act by Av : 

Act = 1/ (JS-2) (57) 

Av = 1 / (KS-4) (58) 

The coordinates of the grid points are given by: 

CTj = (j-2) Act j = 1 ,2,3. . . JS , JS+1 (59) 

v k - (k-4) Av k = 4,5,6. ..KS (60) 

In addition, time is advanced in increments of Ax: 

t. = (i-1) Ax i = 1,2,3...“ (61) 


where the value of Ax is determined by a combination of stability and geo- 
metric constraints which will be discussed later. The upper limit on the 
time counter i is not generally known a priori, but it should be sufficiently 
large that an asymptotic solution is attained. 
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The same values of JS and KS are used in all seven domains. The grid 

columns j=2 and JS correspond to the axial boundaries m=m. and rri. , , . The 

columns j=1 and JS+1 overlap into the adjacent domains and are used to patch 

the solutions together. At the inlet station (which may be either m^ or m^) 

the grid column j=2 is composed of boundary points (to be discussed later) and 

the column j=l is unused. Correspondingly, at the discharge station (either 

m, or m Q ) the grid column j=JS is composed of boundary points and the column 
/ o 

JS+1 is unused. 

The grid rows k=4 and KS correspond to the circumferential boundaries of 
the b 1 ade- to-bl ade passage 0=0^ and 0 u - The grid rows k = 1 ,2 and 3 and 
k = KS+1 , KS+2 and KS+3 are exterior to the computational domain and are re- 
served for the solution in portions of the adjacent blade-to-blade passages 
given by: 


(k-3) Av 

k = 1 ,2 and 3 

(62) 

(k-5) Av 

k = KS+1 , KS+2 and KS+3 

(63) 


Thus, it can be seen in Figure (5) that the r 9 WS k=3 and 4 occupy the same lo- 
cations on the transformed boundary v=0. The rows k = KS and KS+1 correspond- 
ingly occupy the same boundary v— 1 . On the slipstreams these points also oc- 
cupy the same physical locations. This convention was adopted in view of the 
anticipated doub 1 e- va 1 ued solutions which will pertain to each side of the 
blade slipstreams which form the circumferential boundaries of domains 4, 6 
and 7. On the blades, the row k=3 occupies a position corresponding to k=KS, 
and 4 to KS+1, as indicated in Figure (5). 


1 3 

The finite-difference algorithm developed by MacCormack is employed to 
carry out the solution at the interior grid points. This algorithm is based 
on the following second-order approximation to the time derivative 


e . 


i+1 , j ,k 


= e. 


i , j , k 


1 ((. 3 ®.) 

2 U 9x ; 


i+1 , j ,k 


<!f> 


) Ax 


i ,j >k 


where the subscripts i,j,k refer to the discrete values of t, a, and v de- 
fined by Equations (59), ( 60 ) and (6l), and the superscript p indicates a 


( 64 ) 
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FIGURE 5. GRID ROW ORDERING IN VICINITY OF BLADES AND SLIPSTREAMS 





"provisional" (or "predictor”) value. The "provisional" values ai 
in a first or "predictor" step: 


e 


P 

i+1 ,j ,k 


e . 


i >j 


(i£) 

'ax' 


Ax 


» J * k 


Thus Equation (64) can be rewritten as: 


e = 1 ( e + e P . + (— ) 

i+l, j,k 2 v i,j,k e i+l ,j ,k 3 t' 


At) 


i+l , j ,k 


which constitutes the second or "corrector" step. 


It should be noted that the first or "predictor" step given f 
(65) can be carried out consistently with the initial value chara> 
interior point solution, i.e., it only involves the known values » 
and spatial derivatives thereof at time t. as given by Equation (■ 
the second or "corrector" step, Equation (66), requires knowledge 
ary point values at time t. + ^ to complete the evaluation of the si 
ti\>es at the points adjacent to the boundaries at time t. + .. The 
practice, the interior point solution cannot be completed by this 
analogous two-step type algorithms without first completing the b' 
solutions, which will be discussed subsequently. it will be show 
boundary point algorithms used herein can be carried out entirely 
the known solution at t. and therefore, in fact, are executed pri 
terior point algorithm. 


The MacCormack algorithm achieves satisfactory stability wit 
tion of artificial damping terms or numerical filtering procedure 
of first-order, non-centered finite-difference approx i mat i ons to 
derivatives which alternate direction between the first and secor 
resulting solution is. however, considered to be second-order acc 

space and time due to the combination of alternating direction di 

1 3 

in the second or "corrector" step . For example, on the "predie 
spatial derivative may be approximated by: 


(±L) 

ou 


f. . . - f. . 1 , 

1 , J ,k I ,J-1 ,k 


f * J > k 


A a 



Then on the corrector step 


(it) 

v 9 ct ' . 


i+1 ,j ,k 


fP - f p 

i + 1 , j+1 ,k i + 1 , j , k 

Act 


(67b) 


must be used. The same procedure applies to the v derivatives. Since the 
order in which the direction of the differences is evaluated is arbitrary, 
it may be cyclically rotated to avoid imposing a preferential bias in the so- 
lution. The rotation algorithm is illustrated in Table I. The central grid 




- TABLE 
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ROTATION ALGORITHM FOR SPATIAL 




DERIVATIVE EVALUATION 


Time Step 

1 terate 

o Derivative Indices 

v Derivative Indices 

• 

1 

j + 1 > 

j 

k+1 , k 

i 

2 

j > 

J-1 

k, k-1 

i + 1 

1 

j + 1 > 

j 

k, k-1 

i + 1 

2 

j , 

j-1 

k+1 , k 

i+2 

1 

j > 

j-1 

k, k-1 

i+2 

2 

j + l . 

j 

k+1 , k 

i+3 

1 

j , 

j-1 

k+1 , k 

I +3 

2 

j + 1 , 

j 

k, k-1 


point is the j,k point in all cases; thus the circumferential position index 
maintained in evaluation of the a derivative is understood to be k, and the 
streamwise position index maintained in evaluation of the v derivative is 
understood to be j. Note that during any time step, combination of the first 
and second iterates produces a time-split central difference, e.g.: 


2 




i . j 


€> 


p 

) 

i+1 ,j ,k 


f? 


- f p + f _ f 

i+1,j+1,k ?+.1,j,k ? , j , k ?,j-1,k 

2Act 


( 68 ) 
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MacCormack^ has indicated that cyclic rotation should enhance the 
stability of the system. However, in the authors experience it has also been 
found to amplify numerical oscillations in the vicinity of a, shock wave (i.e., 
pressure undershoots and overshoots). The latter phenomenon is not an in- 
stability in the usual numerical sense, but rather ampli.fi cat ion of oscilla- 
tions associated with representation of a discontinuity by a continuous func- 
tion. For this reason, the option to not use the combinations shown in Table 
I for steps i+1 , i+2 and i+3 has been retained in the computer code, at the 
discretion of the user. 


In reference to convergence and stability, the CFL criteria places an 
upper limit on the permissible time step: 


A t < m i n [ {——) 


a (x . - 
- o i+l 


- x . ) 


-). (tit) 


u+a' 


v+a 


a o (y ..' y j3 ) 

( — r • )] 


(69) 


The application of boundary conditions in the case of unequal numbers of blades 
in the two blade rows (which will be discussed later) is based on a phase lag 
in time. Implementation of the phase relations requires a constant value of 
At; therefore an estimate of the maximum anticipated values of u+a and v+a 
must be made to determine the allowable step size At. In addition, the value 
must then be reduced such that the time for one blade of the second row to 
cross a single passage of the first row is an integer number of time steps. 

(The latter constraint obviously does not apply when only a single row is con- 
s i dered . ) 


Inlet and Discharge Boundary Conditions 
and Boundary Point Solution Algorithm 


The present formulation assumes that the periodicity of the flow field 
results entirely from the interaction of the rotor and stator blade rows. Ac- 
cordingly , any nonuniformity, either spatial or timewise, of the properties 
of the flow crossing the inlet or discharge station attributable to inward 
travelling waves or inward convection must be ruled out, since its existence 
would add unsteady components to the flow which are not accounted for by the 
blade row interaction model. Identification of the flow properties which 
propagate by convection and by wave motion is facilitated by recasting the 
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system of equations in character? st i c form. 

The characteristic surfaces formed by the hyperbolic system of partial 
differential equations given by Equations (27) through (30 ) consist of a conoid 
with its base on the x,y plane and within it a stream path which intersects the 
conoid at its vertex. If the vertex is placed at a grid point at time t+At, 
the base covers the domain of dependence of the point at time t. A particularly 
useful approximation to the true characteristic form is obtained by stating 
the system of equations in a reference plane coordinate system which reduces 
the problem to a more tractable two-dimensional (i.e., time and distance) form. 
If the reference plane is normal to both the x,y plane (t=constant) and the in- 
let station (x=constant) and is allowed to translate in the y direction at the 
same velocity as the circumferential component of gas velocity, v, then (as 
will be shown) the system of characteristic lines illustrated in Figure (6) is 
obtained. The lines AO and CO approximate the intersection of the reference 
plane and the true characteristic conoid, and can be interpreted as paths of 
downstream and upstream travelling waves. The line BO is the stream path, 
which also lies in the reference plane. 

It is also pointed out parenthetically that translation of the reference 
plane at the velocity v in effect transforms a rotating frame of reference 
back to an absolute frame. Thus, the inlet and discharge station solutions 
are effectively carried out in an absolute frame regardless of the relative 
motion of the computational domains, and the numerical results are entirely 
independent of this trans format i on . 

Those characteristic lines which originate outside the computational do- 
main at time t and intersect a boundary point at time t+At each represent at 
least one equation which must be replaced by a boundary condition to render a 
determinate solution at the boundary point. Obviously, a subsonic axial Mach 
number, i.e., (u-a) < 0, has been assumed in the construction of line CO in 
Figure (6). Should the axial Mach number at the inlet be supersonic, all 
characteristics would originate outside the computational domain, and accord- 
ingly the complete solution at this boundary could be specified as a boundary 
condition. Correspondingly, a supersonic axial Mach number at the discharge 
station would imply that no characteristics originate outside the computat i ona 1 
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FIGURE 6. CHARACTERISTIC LINES AND GRID POINTS AT INLET AND DISCHARGE STATIONS 






domain and, therefore, no boundary condition can be specified. The present 
analysis is formulated with respect to subsonic axial Mach numbers, but con- 
sideration of supersonic axial flow clearly requires only a minor variation 
of the solution algorithm. In this case, specification of three (3) boundary 
conditions at the inlet will be required (one to replace the wave motion chai — 
acte'ristic AO and two to replace stream path characteristics on BO) . At the 
discharge station one (l) boundary condition will be required to replace the 
wave motion characteristic on CO. Derivation and discussion of the character- 
istic relations and boundary point solution algorithms follows. 

(a) Stream Path Characteristics - The energy equation as given by Equa- 
tion (11) is already in character i st i c form; it can be integrated to yield 
(in terms of the two-dimensional system given by Equations 20 to 2b): 

S = constant on — = — = dt (70) 

Equation (70) states the well known fact that in unsteady flows entropy is 
convected on stream paths. It is evident that the entropy convects inward 
(assuming u > 0) across the inlet station and outward across the discharge 
station. The condition S = constant everywhere upstream of the inlet station 
forms the first boundary condition to be applied at the inlet/ 

The momentum equation, Equation (2), can be rewritten in terms of the gra- 
dients of total enthalpy and entropy rather than pressure: 


— = - VH + V x V x V + TVS 

If VS = 0, as at the inlet, then the curl of Equation (71) gives the vorticity 
transport equation: 

— + Vx (wxV) = 0 


*The value of the entropy at the inlet could, of course, be increased by upstream 
travelling shock waves. However, it is assumed that any shock waves reaching 
the inlet station are sufficiently weak to be considered isentropic. 


2b 
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where the vorticity vector is defined by u = VxV. 

Equation (72) can be expanded and combined with the continuity equation, 
Equation (l), to obtain: 

»(“) = “ (73) 

Dt p p 

In terms of the presently considered two-dimensional system given by Equations 
(20) through (24), the only non-zero component of vorticity is the normal com- 
ponent £ = (3v/8x - 3u/3y) . Accordingly, Equation (73) reduces' to: 

& <f-> ■ 0 

or 

— = constant on — = = dt (75) 

p u v 

Thus, it can be seen that tie ratio of vorticity to density also convects on 

J- 

stream paths. The condition £ = 0 upstream of the inlet station is taken as 
the second boundary condition to be applied at the inlet station. 

To summarize, Equations (70) through (75) apply on the stream path char- 
acteristic, line BO in Figure (6), and have been replaced by the boundary con- 
dition S = constant and C = 0. 

(b) Wave Motion Characteristics - The equation of state. Equation (12) can 
be differentiated with respect to time using the convective operator: 


1 Dp , y Dp j 1 PS 
p Dt p Dt C v Dt 


(76) 


*ln th'e present context this result is restricted to conditions for which 
3S/3y = 0, e.g., at the inlet station. It can also be proven for barotropic 
and constant density fluids. 
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The continuity equation. Equation (27), and the energy equation. Equation (39b), 
can be substituted into the above relation, and the streamwise momentum equation, : 
Equation (28b), then added to and subtracted from the result to obtain: 



where and are identified as the terms appearing on the right hand side of 
Equation (77). Equation (78) applies on the characteristic lines originating 
at points A and C at time t which intersect the boundary point 0 at time t+At 
in Figure (6). The set of equations represented by Equation (78) are commonly 
referred to as the compatibility relations. The member of the set with the + 
sign applies on the line AO at the inlet and will be replaced by an inlet boundary 
condition. The member of the set with the - sign applies on the line CO at the 
discharge station and will also be replaced by a boundary condition. Applica- 
tion of the compatibility relations on line CO at the inlet station and line AO 
at the discharge station to complete the boundary point solution algorithms 
will also be outlined below. 


(c) Modelling of Duct Boundary Conditions - As indicated above in connec- 
tion with the stream path description, it is assumed that S = 0 and £ = 0 at 
the inlet. Th.ese conditions derive from the convective character of Equations 
(70) and (75) and are independent of the physical structure of the inlet duct. 
On the other hand, the propagation of waves across the inlet or discharge sta- 
tion is dependent on the configration of the duct, as ? s we 11 known from 
acoustic theory which ascribes impedance functions to the geometric and ma- 
terial properties of the duct. (Since the present model neglects the radial 
velocity component, radial wave modes and wall impedance properties are coi — 
responding ly assumed to be absent.) 

Two limiting cases have been considered to be descriptive of the duct con- 
figuration at the inlet or discharge station. One case is an Infinite duct, 
i.e., the inlet or discharge station is located in a region of constant cross- 
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sectional area which extends "very far" outward from the inlet or discharge sta- 
tion 1 . Consequently, all outward radiating waves should pass the inlet or dis- 
charge station without reflection in this case, and continue travelling outward 
"forever" on the time scale of the problem. No reflected waves from x -*■ +_ °° 
should ever reach the computational domain. 

The second case represents the opposite limit in which all pressure waves 
are reflected at the inlet or discharge station; this case is termed an open- 
end duct since it corresponds to the type of reflection associated with the 
open-end of an organ pipe. The boundary condition is that the pressure matches 
the plenum pressure outside the duct, namely: 


Ptnlet (y>t) = P -< 


(79) 


^discharge 


(30) 


where the subscripts +°° denote x>x,. , and x<x. , ^ respectively. 

K — discharge inlet 1 


Modelling of the non- re f 1 ect i ve condition for an infinite duct is some- 
what more complex, particularly in regard to the swirling waves produced by a 
rotor-stator interaction. A precise mathematical formulation of the flow field 
solution upstream of the inlet station and downstream of the discharge station 
based on a sma 1 1 -per turbat i on analysis is described in Volume Ml of this report, 
and is herein referred to as the acoustic fai — field model for an infinite duct. 

It accomplishes the desired objective of allowing an arbitrary transient signal 
to radiate outward without reflection and asymptotic attainment of a periodic 
solution with as many harmonic components as can be derived from the 'number of 
grid points spanning the considered boundaries. However, an approximate model 
of the infinite duct condition has also been developed which does not require 
use of the acous t i c far- f ie 1 d analysis.. In the approximate model, the infinite 
duct conditions are derived from the wave-motion characteristics represented 
by Equation (78). Consequently, the discrete acoustic modes are not explicitly 
identified in the approximate model, and their unimpeded transmission across 
the boundary cannot be guaranteed. In the approximate infinite duct model. 
Equation (73) is integrated to yield: 
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+ 0_ 2 ) dt 


constant 


on 


dx _ dy 
u+a v 


dt 


(81) 


where S = constant along the wave path is assumed. Outside the computational 
domain r = constant and b = constant is assumed. The remaining term in the 
integrand of Equation { 8 1 ) , namely a 3v/3y, accounts for the two-dimensionality 
of the actual wave surfaces, as compared to the one-dimensional (helical) sur- 
faces which would result if the swirl component of velocity, v, were constant 
or only a function of x. If 3v/9y is neglected outside the computational domain, 
then the well-known Riemann invariants for the incoming waves are obtained from 
Equation ( 8 1 ) : 

(82) 
( 83 ) 

In this case the subscripts denote the specified values for x +_ °°, i .e. , the 
"ends" of the infinite duct. The two-dimensionality of the outward radiating 
waves at the inlet or discharge station is retained by evaluating the inte- 
grand of Equation ( 8 1 ) numerically: 

t+At 

± u = (~hf ± u ) a + J ("a^ ± ^ dt 

C t 

where: 

x = x - (u+a) At (85) 

3 3 

C C 

Some distortion of the swirling waves at the inlet and discharge bound- 
aries can be expected to result from the above model of infinite duct boundary 
conditions, due to the representation of the incoming waves by the Riemann 
invariants for a one-dimensional unsteady flow. The severity of the distortion 
will depend on the relative strengths of the two wave systems and, therefore, 
should not be serious since the outgoing waves (which produce the swirl) are 
described by a two-dimensional system of equations. In no event should this 


(- 2a Y’ t} + u (y , t ) ) . . - — ^ + u 

' y- 1 inlet y-1 


- u (y - 0) di 


2a 

00 

scharge y-1 
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model generate the standing waves associated with an open duct condition, but 
neither can it be expected to duplicate the perfect wave transmission of the 
acoustic fai — field model described in Volume III of this report. Therefore, 
the approximate infinite duct boundary conditions are considered to offer a 
useful analytical tool for analysis of the aerodynamic performance of intei — 
acting blade rows without Invoking the additional computational complexity of 
the acoustic far-field analysis. 

(d) C ircumferential Velocity Solution - The determination of the swirl 
component of velocity at the boundary points at time t+At can be accomplished 
by use of the circumferential component of Equation (71): 


9_v 
9 1 


9H 9 S 

9y 9y 




' ( 86 ) 


Equation (86) pertains to the flow normal to the reference plane in which the 
wave paths have been identified above. Thus, the swirl component of velocity 
does not exhibit a wave-like behavior in this formulation. Furthermore, Equa- 
tion (86) does not involve any streamwise gradients (if £ is known) and can, 
therefore, be evaluated along the inlet or discharge boundary by the same finite 
difference algorithm employed at the interior points. The swirl component of 
velocity is, of course, implicitly coupled to the axial component and to the 
pressure through the gradient of total enthalpy 9H /9y (even when S = £ = 0) . 


(e) Inlet Solution Algorithm - The inlet solution algorithm consists of 
the stated boundary conditions, namely: S = constant, C =0, and either 

P.™ = constant (open duct) or 2a /(y - l) + u = constant (infinite duct), to- 
gether with the compatibility equation on the upstream wave, line CO in Figure 
(6), given by Equation (78) and the circumferential momentum equation given by 
Equation (86). (This combination of 3 boundary conditions and 2 equations may 
appear redundant since there are only 4 dependent variables, however, the condi- 

, v ' e ! i ■ ' ■ > ■_< : 

t i on C_ co = 0 is a Neumann type boundary condition which only serves to allow 
solution of Equation (86) without knowledge of streamwise derivatives.) 


The boundary condition S_ co = constant is enforced by requiring that two 
thermodynamic properties of the inlet flow which define the entropy be speci- 
fied, e.g., = C v log (p_ co /kp^ co ) = constant. Since the open duct condition 

requires specification of p^ and the infinite duct condition requires a , 
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these two variables have been selected to define S_ =0 and, therefore, k = 
(a^ p ^ ^ j n j 30t f 1 cases. 


— OO ‘ — CO 


The solution algorithms for the open duct and for the infinite duct inlets 
are summarized as follows: 


(1) Open-End Inlet Duct - The values of a_ ro and p_ co are specified, and 
S = £ = 0 is implied. this combination provides the values of p,p and e 


-OO -CO 

at the inlet: 


P = P 


p = yp/a' 


( 87 ) 

( 88 ) 


e = p/ ( (y- Op) 


( 89 ) 


The axial velocity is obtained from the integrated form of the compatibility 
relation, Equation (84), on the upstream travelling wave: 


2 (a-a ) 

u ■ U C + - (F T ) '• + (aQ t + V A 

C 


(90) 


where the characteristic point x , y is located relative to the boundary point 

— C - — 

x,y from: 


x = x - (u -a ) At 

c c c 

y = y - v At 

c c 


and the local sound speed is given by 


(91a) 

(91b) 


a = yp/p 


(92) 


The c i rcumferen t ia 1 velocity component is obtained from f i n i te-d i f fe rence so- 
lution of Equation (86), using S = C =0: 


5v 

at 


3H 

9y 


( 93 ) 


(2) I n f i n i te Duct 1 n i et - The values of a 

ap d S_ oq = £ = 0 is implied. In this case, the internal energy is obtained 


, p and u are specified, 

— OO 1 — OO — OO r ’ 
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from the speed of sound, using the compatibility relation (Equation (84)) 
on the upstream wave: 


e = 


+ T 1 (“..-V (a W c At) j 2 /(y(Y-0) 


2 ' c — - - - c 

The axial velocity component is obtained from Equation (82): 


( 94 ) 
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The pressure and density are given by: 
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— oo 


(95) 
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YP 
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(97) 


The circumferential velocity component is again obtained from Equation (93). 


The pressure is a redundant member of the set of dependent variables, since 
it can always be obtained from the density and internal energy, i.e., p = (y-l)pe; 
however, it is carried in the above presentation of the inlet solution statement 
for clarity. 


The point x , y^ is located by iteration, using linear interpolation to 
determine the flow properties at the point. A change in position of 0.1% is 
used as the convergence criterion. (n principle, the terms multiplied by At 
in Equations (90), (91a), (91b) and (94) should be replaced by the average of 
the values at point C and the new values at point 0; however, examination of 
the solution a posteriori indicates that the numerical error involved is too 
small to warrant the additional computational complexity. 

(f) Discharge Station Solution Algorithm - The considerations pertaining 
to formulation and statement of the discharge station boundary conditions and 
solution algorithms are similar to those outlined above. However, as indicated 
in Figure (6), points A and B originate within the computational domain, and 
point C falls outside. Therefore, the compatibility relation pertaining to 
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point C is replaced by a boundary condition given by either Equation (80) or 
( 83 ). Equation (70) is evaluated on the stream path B in this case, using 
linear interpolation. The discharge station boundary conditions and solution 
algorithm are summarized by: '■ ■ 

(1) Open-End Discharge Duct - Value of p^ is specified. 


P = p b (p/p b ) ,/Y 
e = p/((y-l)p) 
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(2) Infinite Discharge Duct - Values of a and u are specified. 
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In both cases, the c i rcumferent ia 1 velocity component is given by 


(104) 
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and the characteristic point x , y is located with respect to the discharge 
boundary point x,y from: 


3 v 
3t 


3H 

3y 


3S 

3y 


- u ? 


x = x - (u + a ) At 
a ' a a 


( 107 ) 


32 



y 


v At 
a 


( 108 ) 


Equations (99) and (10*0 utilize Equation (70), i.e., S = S^. Solution of 
Equation ( 1 06) is carried out using a first order upstream difference (i.e., 
j, j-1) to evaluate the a derivatives necessary to compute £ at the boundary 
points . 


Blade Surface and Slipstream Boundary 
Conditions and Solution Algorithm 

(a) Blade Surface Boundary Conditions - The boundary condition at the 
blade surface is simply impermeability of the surface, which requires that the 
component of velocity normal to the surface vanish. The blades are assumed to 
be thin and have sharp leading and trailing edges, as is typical of high speed 
compressor and fan blades. Therefore, the Kutta condition which requires the 
pressure to be continuous and finite at the trailing edge, is applicable. This 
also implies that the slipstream (i.e., vortex sheet) which emanates from the 
blade must originate at the trailing edge. Since the blade leading edge is 
sharp, and the incidence angles are not expected to be large, the streamline 
which wets the blade surface is assumed to intersect the leading edge. Accord- 
ingly, the leading edge pressure is also required to be finite, but not continuous. 

(b) Slipstream Boundary Conditions - The boundary conditions pertaining 
to the slipstream are similar to those for the blade surface in that the slip- 
stream is impermeable, but are dissimilar in that the slipstream is non-rigid. 

It can be shown from application of the conservation form of the governing equa- 
tions at an impermeable contact surface that the pressure must be continuous 
across the slipstream, and that the component of velocity normal to the slip- 
stream surface must also be continuous and equal to the surface velocity. 

(Thus, in a frame of reference moving with the surface the normal component of 
velocity must vanish at the surface.) However, the tangential component of 
velocity may be discontinuous across the slipstream, as well as the density or 
other thermodynamic properties. These jumps (discontinuities) in flow proper- 
ties result from unsteady variations in the work performed by the blades, or by 
differences in shock-produced losses on either side of the blade under steady 
conditions, for example. It is emphasized that the conservation form of the 
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equations admit the existence of these jumps across a slipstream, but their 
magnitudes are not derivable from appl icat ion of the conservation form of the 
governing equations to a surface of discontinuity. By contrast, the magnitude 
of the jumps across a shock wave derive from the Rankine-Hugoniot relations. 
Therefore, the slipstreams cannot be expected to evolve from a finite difference 
solution of the governing equations in the same way that shock surfaces are 
"captured”. Consequently, if the slipstream jumps and the corresponding slip- 
stream motion are to be resolved accurately, the slipstreams must be explicitly 
recognized as surfaces of discontinuity, as they are in the present formulation. 


(c) Restatement of Governing Equations - The solution algorithms used at 
the blade surface and slipstream points are closely related and, therefore, will be 
derived for the more general case of the slipstream points. The result will then 
be specialized to the blade surface points and the leading and trailing edge points. 


A surface-oriented coordinate system (x,y,t), as sketched in Figure (7), 
is defined by: 
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where x is 
mal to the 
ponents ; n 


the curvilinear distance along the surface, y is the distance nor-* 
surface, and <j> is the angle between x and x. The velocity corrt- 
this system (u,v) are correspondingly defined by: 


u = ucostj) - , vs I n<j) 

v = us i n<f> . + vcos<j) , , 

In this system. Equations (27), (28b), (29b) and (30b) become: 
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PROJECTION OF SLIPSTREAM AT TIME 
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FIGURE 7. CHARACTERISTIC CONSTRUCTION AT A SLIPSTREAM POINT 
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where 
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Combination of Equations (113) and (116) and the equation of state yields: 
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Equations (11^) and ( 1 1 8 ) can then be combined to obtain a pair of compati- 
bility relations analogous to those previously discussed in connection with 
the inlet and discharge boundaries: 
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The variables Q^and Q 2 are defined by identifying the right hand side of Equa- 
tion (120) term by term with the right hand side of Equation (119). 


Implementation of the numerical solution of Equation (120) is facilitated 
by introduction of an additional coordinate t ransformat ion from the curvi- 
linear (x,y) system to a series of local Cartesian systems (S,n) each of which 
is tangent to a grid point on the slipstream (or blade surface). The (C,n) 





system is shown schematically in Figure (7). Note that in the (£,n) system, 

R”^ = 3<j>/3t = 0. The velocity components (u,v) at the grid point at which 
the local (£,n) system is defined are unchanged by this transformation; how- 
ever, their values at adjacent grid points must be evaluated with respect to 
the angle <j> at the subject grid point. The overall effect of these transforma- 
tions is to make evaluation of Equation (120) closely approximate impingement 
of a one-d i mens i on-a 1 acoustic wave on a surface which is moving at a velocity 
v. Although all twO-d i rnens i onal terms are in fact retained, they can be viewed 
as "corrections" to a more familiar one-dimensional solution. 


The boundary point algorithm is completed by stating the streamwise momentum 
equation and the energy equation in the curvilinear (x,y) system which follows 
a stream path along the moving surface. The magnitude of the velocity vector 
along the surface is: 
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The streamwise momentum equation is obtained by multiplying Equations (114) 
by v and (115) by u and adding: 
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(123) 


Equation ( 1 1 6 ) integrates exactly to: 


c ' ds _ dx _ dy . 

S = constant on — = — = = dt 


(124) 


Equation (122) is not an exact integral, and, in contrast to the compatibility 
relations, i.e., Equation (120), the integrand on the right hand side of Equa- 
tion (122) includes a leading term (the pressure gradient) which cannot be 
adequately approximated by its initial value on the path of integration. Note, 
however, that Equation (30a) represents a combination of the momentum and energy 
equations; therefore since Equation (30b) has been used to represent conserva- 
tion of energy. Equation (30a) can be used to represent conservation of momentum 
in lieu of Equation (122). Either form should be equivalent, but it will be 
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seen below that Equation (30a) offers an advantageous form for numerical 
evaluation. Equations (27) and (30a) can be combined to yield: 


8H 3H_ 3H_ J_ 3p 

3t U 3x V 3y p 3t 

In the present coordinate system this becomes 
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where the convective derivative is defined above by Equation (123). In most 
applications for unsteady flows, Equation (126) is awkward to evaluate because 
the time derivative of pressure is not known a priori. However, in the pre- 
sent formulation, the pressure at time t^ (point 0 in Figure 7) is determined 
from the solution of Equation (120), which is uncoupled from Equation (126) 
and, therefore, can be evaluated prior to (126). Thus, the term in question 
can be accurately approximated by: 


1 iP 
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Accordingly, Equation (726) integrates to: 
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( 128 ) 


The constants indicated in Equations (124) and (128) are determined by 
evaluating S and H , respectively, at a distance As = q At upstream of the point 
Q. (in Figure 7). Although neither Equations (126) or (122) are exact integrals, 
as is Equation (124), the time derivative of pressure along line 00 (in Figure 
7) can be more accurately represented than its spatial derivative at point 
B or E (which must be interpolated). Numerical experimentation has shown 
that integration of Equation (126) yields correspondingly more accurate re- 
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suits than integration of Equation (122). 


(d) Blade Surface Solution Algorithm - Equations (120), (124) and (128) 
constitute the system of governing equations as applied on the blade surface 
and slipstream points. Therefore, with reference to the geometry indicated in 
Figure (7), and reverting to the subscripts SL and u to denote lower and upper 
boundary surfaces, respectively, the solution at a blade surface point 0 is given 
by: 
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Determination of the total internal energy follows from the total enthalpy, 
pressure and density by definition, c.f. Equation (26). 


(e) Slipstream Solution Algorithm - The solution at a slipstream point 
is given by: 
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The overall procedure for imposing the boundary conditions along the blade 
surface and slipstream points is shown schematically in Figure (8), The time 
axis projects vertically out of the page in this figure. The dashed lines re 
present the intersections of the reference planes and of the stream path with 
the axisymmetric stream surface (i.e., the x,y plane) during a time step At 
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(136d) 

(137a) 

(137b) 
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FIGURE 8. COMPARISON OF BLADE AND SLIPSTREAM CHARACTERISTIC SYSTEM 



(neglecting the motion of the slipstream during the interval At). The values 
obtained for v on the slipstream are used to locate the new position of the 
slipstream for the next time step as follows: 

x = x (vsincj>) At (138) 

no 

y = y + (vcoscf>) At 039) 

no 

where (x , y (x ,t)) are the coordinates of the slipstream point at which 
000 

the above described solution is obtained, and (x , y ) are its new coordinates. 

n n 

Linear i n te rpo 1 at i on is then used to find the new y coordinate, y (x , t+At) 
r 00 

of the intersection of the slipstream and the grid column located at x=x , which 

o 

gives the new coordinates of the slipstream grid point at the next time step. 


(f) Trailing Edge and Leading Edge Points - At the trailing edge, the 
system given by Equations (129) through (132) is applicable subject to the 
constraint that p =p . This condition is satisfied by iterating the angle <{> 
at the trailing edge, which, unlike the general blade surface and slipstream 
points, is not known a priori. The blade trailing edge is assumed to be sharp, 
but not necessarily cusped; the angle <|> is taken as a weighted combination of 
the blade surface angle and the slipstream angle, as shown schematically at 
points A and B in Figure (9a). 
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(Values of k-^ in the range y <_ k.^ <_ -j have been found to provide satisfactory 
accuracy and stability of the trailing edge point solutions.) Therefore, the 
trailing edge point solution is carried out at points A and B like any other 
blade surface point solution, except that the slipstream angle is deflected 
until equal pressures (to within an error tolerance of 0.0001% of the mean 
pressure) are obtained on each side of the surface. The trailing edge and the 
slipstream emanating from It thus can be considered as together forming a con- 
tinuous deformable surface. Note, however, that one value of <j>_ E will be ob- 
tained on the upper surface of the blade and a second value on the lower sur- 
face, if the trailing edge itself is not actually cusped. Therefore, the direc- 
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tion of the velocity vector at the trailing edge is multiple-valued. Spe- 
cifically, 5 values of the flow angle are i dent i f ied, at 6 points having a 
common location; the trailing edge surface angles at points E and F, the 
mean angles at points A and B, and the (single) slipstream angle at points 
C and D. Since the Kutta condition requires the pressure to be continuous 
at the trailing edge, the change in pressure associated with the angular dif- 
ferences between these 6 points is neglected. Therefore, the solution is 
carried out at point A as indicated above, and then solutions are defined 
at points D and E which only differ from that obtained at point A in the 
direction assigned to the velocity vector. (i.e., the pressure, density and 
magnitude of the velocity vector at points A, D and E are equal. The direc- 
tions of the velocity vectors are parallel to the angle of each of these 3 
points.) The same procedure is used to carry out the solutions at points B, 

C and F. Note that velocity vectors obtained at points C and D are, there- 
fore, parallel but can differ in magnitude. The solutions at points A and B 
are used in the computations at interior points directly above and below the 
trailing edge, whereas the solutions at C and D are used to determine the 
slipstream point solutions downstream of the trailing edge, and those at 
points E and F are used to carry out solutions at the blade surface points 
upstream of the trailing edge. 

The leading edge configuration is sketched in Figure (9b). In view of 
the singular character of a sharp leading edge, a somewhat more complex solu- 
tion procedure than employed at the trailing edge is necessary. A total of 
five flow angles, and a corresponding number of solutions, are again identi- 
fied at the leading edge. With reference to Figure (9b), the points A, B, C, 

D and E are all located at the same physical position, however, the flow angle 
at point A is the angle of the stream path intersecting the leading edge, and 
the flow angles at points B and C are a weighted combination of the surface 
angles of the upper and lower surfaces, respectively, and the stream path 
angle : 


^LE ^LE ^surface + ^ ^LE^ ^stream path 

(A value of k^_ = 1.0, which makes the flow angles at points B and C tangent 
to the actual blade surface angles, has been found to be satisfactory for very 
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slender blades, but k L£ <l .0 may be appropriate for other cases). The flow angles 

at points D and E are the mean of <f> , and * , 

LE stream path. 

The leading edge solution procedure is begun by varying the stream path angle 
until the pressures at point A obtained from Equation (130a) and (130b) are 
matched (to within an error tolerance of 0.0001% of the mean pressure). As part 
of this step, the instantaneous locations of the origins of the stream path and 
the wave characteristics intersecting the leading edge are determined. The solu- 
tion at point A is then completed using Equations (129), (131) and (132) in the 
same fashion as employed at any blade surface point. Note that a s i ng le-val ued 
solution will be obtained at point A. Next, pressures at points B and C are 
obtained from Equations (130a) and (130b) using the specified angles at these 
points (from Equation 141). Since points B and C are, in fact, coincident with 

JL 

point A, the solutions must have the same entropies” at any instant. The total en- 
thalpies are determined from Equations (132a) and (132c). Therefore, the 
densities and velocities at these points are found from the equation of state, 
the definition of total enthalpy and the specified angles. Finally, solutions 
at points D and E are obtained by averaging the pressure and magnitudes of the 
velocity vectors at points A and B and at points A and C, respectively. The 
directions of the velocity vectors are also averaged. The densities at these 
points are then calculated from the entropy at point A. 

Summarizing, at the leading edge a single value of the entropy and total enthalpy 
is obtained, but five values of the pressure, density and velocity vector are 
determined. The solution obtained at point A in Figure (9b) is subsequently 
used in computations at the interior grid point upstream of point A. The 
solutions at points B and C are used for the adjacent blade surface points on 
the upper and lower surfaces, respectively. The solution at point D is used in 
computations at the interior point above the leading edge, and that at point E 
for the interior point below the leading edge. 


*The possible existence of a shock wave at this point due to the flow deflection 
is neglected in this connection. The oblique shock entropy increase could, of 
course, be computed, if appropriate, by suitable modification of the code. 
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( g ) Intersection of Slipstream and Dis charge Boundary - Finally, s pec i a 1 
consideration must be given to the grid points at the intersections of the slip- 
streams and the discharge boundary, since these grid points lie on two boundaries 
of the computational domain. It should be emphasized in regard to these points 
that all the characteristic points (A, B, D and E in Figure 7) which influence 
a slipstream point (0 in Figure 7) 1 i e upstream of the subject point (assuming 
u > 0) . The numerical domain of dependence of a slipstream point, therefore, 
extends to the adjacent grid points only through the linear interpolations 
necessary to evaluate variables and derivatives at the characteristic points. 

This dependence represents the only mechanism by which those slipstream points 
which lie on the discharge boundary are affected by the discharge boundary con- 
ditions that are explicitly enforced at all other discharge boundary points. 

In connection with these same special slipstream points, it should be 
pointed out that if the slipstream angle <p at the discharge boundary becomes 
suff ici ently- large (i.e., the included angle between the slipstream and dis- 
charge boundary is sufficiently acute) it is possible for the characteristic 
point to be located outside the computational domain. The physical meaning 
of this occurrence is that the slipstream solution explicitly depends on data 
downstream of the discharge boundary and, therefore, the boundaries of the 
computational domain do not encompass the numerical domain of dependence of 
the solution. Theoretically, this condition represents a limitation on the 
applicability of the present formulation. In practice, it represents an ex- 
treme condition which has thus far only been encountered during transient 
phases of a solution; use of linear extrapolation of data from the set of 
grid points adjacent to the discharge boundary to the character i st i c point has 
proven successful in these instances, although it temporarily violates the CFL 
convergence criterion. 


Periodicity Condition 

In the analysis of an isolated circular (or infinite) cascade of blades 
in a uniform free stream, it is clear that the solution for each blade-to- 
blade passage will be identical. The solution for the complete cascade in this 
case will be steady (in the frame of reference of the blades) and have an 
angular period of 2ir/N. Enforcement of this periodicity condition is straight- 
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forward; the solution on the exterior grid rows k=1 , 2 and 3 (see Figure 5) 
can be equated to those on interior grid rows KS-2, KS-l and KS, and those on 
exterior rows k=KS+1 , KS+2 and KS+3 can be equated to those on interior rows 
k=4 , 5 and 6. It may be noted that in domains 1 and 2, grid rows 1 and 3 
and KS+1 and KS+3 are superfluous, since the lateral boundaries of these do- 
mains are in reality composed of ordinary interior points. However, the so- 
lutions on grid rows 1, 2 and 3 and KS+1, KS+2 and KS+3 are required to carry 
out the slipstream point solutions in domains 4, 6 and 7 as well as the blade 
leading and trailing edge point solutions. (Compare Figures 5 and 8.) 

In the more general case of a stage composed of interacting blade rows 
having unequal numbers of blades, the fundamental angular period of the com- 
plete cascade solution is 2 tt/AN, where AN is the difference in the number of 
blades. Furthermore, the flow pattern rotates with an angular velocity which 
is, in general, a multiple of the wheel speed. ^ (In the special case of an 
equal number of blades the flow pattern rotates at "infinite" speed and the 
angular period 2-?t/N is recovered.) Numerical representation of the periodicity 
condition pertaining to the conventional configuration, consisting of a pair 
of blade rows with the larger number of blades in the second row, has been ac- 
complished by formulation of a cyclic procedure for equating the solution on 
the exterior grid rows identified above (i.e., k=1 , 2, 3 and KS+1, KS+2, KS+3) 
to that on corresponding interior rows (i.e., k=KS-2, KS-1, KS and 4, 5, 6) at 
an earlier time through a set of appropriate phase relations. The same pro- 
cedure is applied on the grid columns at the interface between domains 4 and 5. 

An illustration of the nature of the cyclic procedure devised to enforce 
the periodicity of the solution can be accomplished through use of the following 
simplified configuration. Consider first a stage having three rotor blades and 
three stator blades. For the present discussion, let the stator precede the 
rotor. This con f i gu rat i on is shown in Figure (10a) in both axial and cascade 
projections. At time t all rotor and stator blades are aligned, whereas at 
time t +At the rotor has moved through a fraction of a revolution, and none 
of the blades are now aligned. It is clear in this case that the geometric 
boundary conditions (which determine the angular periodicity of the flow through 
the stage) are identical in each b 1 ade-to-b 1 ade passage at any time. Accordingl 
the flow field in each passage should be identical, since none of the boundary 
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FIGURE 10a. ILLUSTRATION OF CYCLIC ALGORITHM FOR STAGE WITH EQUAL NUMBER OF 
BLADES IN STATOR AND ROTOR (N. = 3, No = 3) 



conditions distinguish one passage from the next. In this case the solution 
along an exterior grid row 6 can be equated to that along the interior grid 
line 3 and similarly that along exterior line a can be equated to that along 
interior line y, at any instant. Consider now the case with three blades in 
the stator and four blades in the rotor as shown in Figure (10b). At time t 
rotor blade 2 is aligned with stator blade b, whereas at time t +At rotor blade 
3 is aligned with blade c. In this case the geometric conditions pertaining to 
the passage between blades a and b are obviously different from those for the 
passage between blades b and c at any time. However, it may be noted that ' 
those pertaining to passage be at t +At are precisely the same as those which 
pertained to passage ab at the previous time t Q . Therefore, the flow condi- 
tions along exterior grid line <5 at time t +At can be equated to those along 
interior grid line $ at the earlier time t . However, in this case those along 
exterior grid line a at time t +At cannot be equated to those occurring in pas- 
sage ab at time t Q , but must be equated to those occurring along line y at an 
earlier time. Thus, a phase shift is introduced in application of the lateral 
boundary conditions. 


A similar procedure is used to define boundary values along the interface 
between domains 4 and 5. However, in this connection it is pointed out that 
sufficient data must be stored along this interface to provide information for 
a maximum period corresponding to the blade passing frequency of the first row 
(i.e., the row with the smaller number of blades). During this period the 
relative angular positions of the two domains will shift by 2 tt/N ^ . In addition, 
domain A wi 1 1 itself span an arc of 2tt/N^ ; therefore data covering a total arc 
of 2 (2 tt/N^ ) must be available. The boundary data is stored for one blade-to- 
blade passage on either side of the central passage which forms the computational 
domain, i.e., a total of three passages. Thus, for domain 5 the stored data 
spans the arc 3(2 tt/N 2 ) . An upper limit on the ratio of number of blades re- 
sults, namely: 3/N 2 >_ 2/N^ is required. The permissible number of blades in 

the second row is, therefore, bounded by: 



The high speed fan configurations to be considered later in this report have 
stator to rotor ratios of the order of 1.05 and 1.12, which are within these 
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FIGURE 10b. ILLUSTRATION OF CYCLIC ALGORITHM FOR STAGE WITH UNEQUAL 
BLADES IN STATOR AND ROTOR (N. = 3, N,*4) 



limits. Should the upper limit prove restrictive for other configurations, the 
stored boundary data for the interface of domain 5 with domain A could be ex- 
tended to span an arc 5(2 tt/N 2 ), for example, by including two passages on 
either side of the central passage, which would increase the upper limit to 
N 2 /N.| 5/2. Some modification of the computer program described in Volume 

It of this report would, however, be required. 

, A more detailed exposition of the development and application of the phased 
boundary conditions is contained in Appendix I of this report. In particular, 
it is demonstrated that the required boundary information from adjacent blade- 
to-blade passages is generated as one passage of the first row crosses one 
blade of the second row, and that the desired periodicity is attained asymp- 
totically in time. 


Initial Cond i t i ons 

The initial conditions need be considered only to the extent that they bear 
on the asymptotic limit in time, since only the asymptotic solution is of 
interest. Several observations relative to this point should be emphasized: 

(a) The initial data is necessarily approximate, at best, since determina- 
tion of the exact solution is the objective of the calculation. The initial 
transient solution is associated with the difference between the initial data 
and the exact solution, and consequently the time needed to attain an asymp- 
totic solution may be expected to diminish as the accuracy of the initial data 

i s i mproved . 

(b) The initial data can be approximate in the sense that it does not 
satisfy the boundary conditions, or it does not satisfy the governing equa- 
tions, or both. If the initial data does not satisfy the governing equations, 
the resulting transient solution has no physical relevance; only the asymptotic 
limit is meaningful. On the other hand, if the initial data satisfies the 
governing equations but does not satisfy the imposed boundary conditions (i.e., 
it satisfies some other set of boundary conditions) then the resulting tran- 
sient solution is physically relevant; it represents the response of the system 
to an impulsive change in boundary conditions. The computer program has not 
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been structured to handle such impulsive changes in boundary conditions, but 
the ability of the formulation to do so is simply pointed out. 

(c) As a consequence of the present formulation of inlet and discharge 
boundary conditions, existence of an asymptotic solution for arbitrary initial 
data or for arbitrary boundary conditions cannot be guaranteed. Obviously, 

it is possible to specify inconsistent boundary conditions, such as an unattain 
ably high pressure ratio across a stage, or initial data that generate a tran- 
sient solution which violates certain underlying assumptions, such as u > 0 
at the inlet and discharge station. An asymptotic solution cannot be attained 
in such cases (i.e., the computer program "bombs"). However, this problem 
should not be any more disturbing than the failure to reach a stable operating 
point in an actual turbomachine due to surge, for example. In the present 
formulation, possible non-existence of an asymptotic solution is the penalty 
incurred by modelling as closely as possible the actual wave mechanics of the 
inlet and discharge flows without modelling the entire starting process by 
which a stable operating point is reached. On the other hand, if existence of 
an asymptotic solution is demonstrated for a particular set of boundary condi- 
tions and initial data, then uniqueness of the numerical solution necessarily 
follows. If it were non-unique, the numerical solution would drift through 
an endless succession of states, since each time step is a new initial value 
problem with a perturbation of the data provided by the round-off error. In 

other words, convergence of the solution guarantees its uniqueness. 

(d) It should be apparent from the preceding discussion that the initial 
data can bear on the existence of an asymptotic solution, but the extent to 
which an asymptotic solution depends on the initial data is more difficult to 
define. It is conjectured here, on the basis of experience thus far, that 
variations in the initial data within those bounds for which a solution exists 
have no effect on the asymptotic solution. It is clear that the solution at 
any instant during the transient oscillation through which the flow proceeds 
can be regarded as an initial condition leading to the same asymptotic solu- 
tion. Although this is not a conclusive proof of the conjecture, it tends 

to support the limited observations on which the conjecture is based. In 
connection with this same question, it should be recognized that application 
of the inlet and discharge boundary conditions formulated herein to a duct 
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flow without, any turbomachinery would not be sufficient to define a unique so- 
lution, as is frequently the case in inviscid flow problems. An infinity of 
solutions, ranging from no flow to choked flow, could satisfy the imposed 
boundary conditions, and only the initial data would determine which solution 
(if any) would be obtained. It is the addition of the rotating blade row and 
the application of a Kutta condition to each airfoil which provides the addi- 
tional constraint that defines a unique solution. From this viewpoint, the 
role of initial conditions should be irrelevant, as long as a solution can be 
obta i ned . 

In the present formulation, initial conditions can be specified in either 
of twQ ways. If no previous information is available. (a "cold start"), the ini 
tial data for the entire computational domain is approximated from the inlet 
and discharge boundary conditions, (supplemented by u_ oo in the open duct case) 
and an initial value of the swirl angle, tan ^ v/u, at the inlet. In this 
case the inlet pressure and entropy are imposed in the domains preceding the 
rotor. The density is computed from the equation of state and the meridional 
velocity from the inlet mass flow rate (pu 2 Trrb) _ ro . The circumferential ve- 
locity component is approximated by maintaining the inlet value of the tangent 
of the swirl angle, (v/u), throughout these domains. Within the rotor, a 
linear increase in pressure from the inlet value at the leading edge station 
to the discharge value at the trailing edge station is assumed. An entropy 
gradient may also be imposed across the rotor corresponding to anticipated 
shock losses. Downstream of the rotor, the pressure and entropy obtained at 
the rotor trailing edge station are assumed to prevail. The density and 
meridional velocity are again obtained from the equations of state and inlet 
mass flow rate, respectively. Downstream of either blade row, the slipstream 
is assumed to 1 i e on a continuation of the trailing edge camber line, and the 
circumferential velocity component is approximated from the slipstream angle 
rather than the inlet swirl angle. The blade surface and s l i pstream .point 
routines are then activated for one time step to satisfy the boundary condi- 
tions on these surfaces. 

If a previous solution, employing the same grid structure, is available 
(i.e., a "restart") it may be used as initial data, with either the same or 
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revised boundary conditions. However care must be exercised in this regard 
when carrying out a computation using the phased boundary conditions. The 
elapsed time must always be measured from the "cold start" since the cyclic 
procedure uses stored boundary data which is identified by a time counter that 
cannot be altered when "restarting" the computation. 

Finally, it is pointed out that an initial transient solution which is 
physically irrelevant can excite slipstream oscillations which are sufficiently 
violent to abort the computation. Therefore, provision has been included to 
utilize a "small disturbance" type slipstream approximation during the initial 
transient phase. In this case, the slipstream point algorithm is carried out 
in its entirety, but the resulting solution is applied on the original s 1 i p— 
stream contour, which is held fixed. The computation of the slipstream motion 
can be restored after the initial transient phase of the interior solution has 
decayed. This option can be regarded as a "two phase" initialization procedure 

BOUNDARY LAYER AND WAKE ANALYSIS 
Motivation and Approach 

The requirement for evaluation of viscous effects in compressors and/or 
fans arises from a number of sources: (a) the displacement effect of the 

boundary layers on the blades, and on the hub and casing walls constricts the 
inviscid flow area and thereby alters the inviscid solution, (b) the viscous 
drag is largely responsible for the loss of total pressure through the stage 
(in a shock-free flow it is solely responsible), (c) passage of one blade row 
through the wakes of a preceding row contributes significantly to the unsteady 
lift of the second row and the associated acoustic properties of the stage, 
and (d) in the extreme, the inviscid portion of the flow is entirely engulfed 
by turbulent eddies. In the last case, the physical and mathematical character 
of the flow is fundamentally changed, and the analytical model must be reformu- 
lated accordingly. In the present approach, as described in the preceding sec- 
tion, it Is assumed that the flow is predominantly inviscid, and viscous ef- 

*Details of these procedures are contained in Volume II of this report. Note 
that the stored data is the non-dimensional form defined by Equation (^9) . 
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fects can be modelled by standard boundary layer theory. Thus, the situations 
in which viscous effects predominate (i.e., separation or fully turbulent flow) 
are excluded. However, approximate representation of the first three effects 
noted above is incorporated in the present model. 

The displacement interaction between boundary layer and inviscid flow under 
steady conditions has been, and remains, a subject of extensive research. Prac- 
tical computational methods which are applicable when separation (i.e., flow re- 
versal within the boundary layer) occurs have not been firmly established as yet 
A less satisfactory situation exists in an unsteady interaction, for which even 
the definition of separation is not universally agreed upon. Unfortunately, the 
pressure gradients in compressors are usually adverse and, therefore, separation 
is a significant problem. 

The intent of the present effort is to provide only an approximation to 
the boundary layer displacement effect on the blade surfaces; therefore re- 
course is made to standard steady boundary layer representations based on zero 
pressure gradient and local similarity concepts with heuristic corrections to 
account for regions of separated flow. If separation is indicated on this 
basis, a more detailed analysis of the boundary layer is clearly required. 


Quas i -Steady Approximations 


The use of steady boundary layer representations implies that either the an- 
ticipated time scale for variations in the inviscid flow is much larger than 
the characteristic boundary layer response time, or that the anticipated un- 
steady components of the viscous flow solution are small compared to the 
steady components. The first condition can be examined by comparing the time 
for a diffusion wave to traverse the thickness of the boundary layer with the 
period between rotor blade passings. The speed of a diffusion wave is of the 

order y/(p6) and hence the time to traverse the boundary layer is of the order 
2 

6 p/p. Therefore, the first condition requires: 
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The maximum boundary layer thickness is on the order of 


6 _ 

c 


- 1/2 

5 R lami nar 

— R turbulent 


(iMO 


where R = puc/u is a Reynolds number based on chord length and characteristic 
(mean) values of velocity, density and viscosity in the inviscid flow. There 
fore. Equation (1^3) becomes: 


(2^) (_iL) 



25 laminar 


turbulent 

y 


045 ) 


(The laminar value is only included for the sake of comparison; a turbulent 
boundary layer is expected to prevail over most of the blade under typical 
operating conditions.) 


Typical values of the parameters appearing in Equation (145) are of the 
order of 2tt/N - 1/10, u/ftr - 1 and r/c - 10. Therefore, the required inequality 
is not generally satisfied; i.e., the response time of the boundary layer is 
longer than the period between blade passings. Consequently, justification of 
the assumption of a quasi-steady viscous solution must be predicated on exis- 
tence of relatively small amplitude unsteady disturbances in the inviscid flow. 
Thus, it is consistent with use of acoustic or linearized inviscid flow theory, 
but not necessarily consistent with the present use of nonlinear theory. The 
possible requirement for an unsteady viscous flow analysis to correctly repre- 
sent the vi scous- i nvi sci d interactions in rotating turbomachinery should be 
recogn i zed . 


Blade Boundary Layer Displacement 
Thickness and Shear Stress 

(a) Basic Equations - In the present context, the principal interest in 
the boundary layer solution resides in determination of the displacement effect 
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on the blade surfaces and the initial conditions for the wakes which emanate 
from the trailing edges. Therefore, the well known momentum integral equa- 
tion (cf. , Reference 15, Chapter VIII) is adequate to provide the required 
information: 

[20 + 6*3 (146) 

where x is the distance along the blade surface measured from the leading edge, 

u refers to the magnitude of the velocity vector in the inviscid flow at the 
e 

blade surface (in the frame of reference of the blade row), p is the corre- 
sponding density, and density variations within the boundary layer as well as 
the effects of variable stream sheet thickness and radius are neglected. The 
shear stress, displacement thickness and momentum thickness are defined by: 
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( 149 ) 


where u is the velocity within the boundary layer. As will be shown below 
Equation ( 1 46 ) can be transformed to the fol lowi ng form for both laminar and 
turbulent boundary layers: . 


^ = ar- b 

dx 


( 150 ) 


where the parameters A and B are approximate functions of the shape factor, 

JL 

6 /0, pressure gradient, Reynold's number , etc. , which will be derived below. 
Numerical integration of Equation (150) by the simple Euler formula provides 
the displacement thickness at the grid points spanning the blade surface: 
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6* +1 = <5* + ( AR " B )j / 051) 

Integration is begun midway between the leading edge and the first surface 
point to avoid the singularity at the leading edge (x = 0) . Evaluation of 
the shear stress is outlined below. 

(b) Laminar Flow - The solutions of the Falkner-Skan equation^, i.e., 

.2 

f + ff + B(1-f ) = 0 (152) 

17 

derived by Hartree (cf. Reference 18 for a complete derivation and discus- 
sion of this equation) form the basis of the laminar flow approximations which 
are utilized up to the transition point. In this case 


f = f(n) 

( 153 ) 



/m+1 U ei 

n i= y (- 5 — p — ) 

(15 1 !) 

yx 


r - 2m 
3 m+1 

(155) 

u = u f (n) 
e 

(156) 


Accordingly, 

i 

T 0 = u e p ~r) f (o) 057) 

yx 

and, from Equation (146) with du /dx = 0: 

e 

f" (o) R _i (158) 

dx J 

Therefore, by comparison with Equation ( 1 50 ) , the values of A and B for a 
laminar flow are identified as 

A = V (2 r->* f " <°) 0 59) 

B = 1/2 (160) 
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curve-fit as a function of the pressure gradient parameter 3 for the attached 
flow branch, 3 ~ .1988: 

JL 

= 2.16 + 1.86 exp(-7.367 (3 + 0.1988)) (161) 

f ”*"* ( o ) = 0.79 (3 + .1988)* 525 (162) 

•* ' JL 

Separated flow solutions are obtained for 3 < - .1988. For this case 6 
and f (o) < 0. The condition 

JL 

-- f~ (o) = 0.0 063) 

is used to permit continuation of the integration of Equation (150) through 
the separation region. However, the result should be viewed with caution if 
separation is indicated, as a more detailed viscous solution is required for 
this case. 


(c) Transition - The criterion for transition from laminar to turbulent 

1 8 

flow is based on a study by Pretsch described in Reference (15), Chapter XVII. 
Transition is considered to occur when: 


R rs* > R 

J critical 
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where R ~ puS /\x. The critical Reynolds number is defined by: 


60 


for 3 < - .1988 


R g * = l 60 + 600 (1 + — |gg ) for -.1988 £ 3 < 0 (165) 

critical * J 


660 + 120003 


for 3 > 0 



(d) Turbulent Flow - A power-law velocity profile is assumed to be de- 
scriptive of the turbulent boundary layer (cf. Reference 15): 
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In this case 
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6 ‘ (l+n) S (2+n) (,68) 
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This type profile was experimentally derived for pipe flows by Nikuradse , 
for which the corresponding shear stress law was proposed by Schlicting 
(Reference 15, Chapter XX): 


x = pu 2 .0225 R* 
o e 6 


(169) 


where R^ = pu6/p. Substitution of these results into Equation (lk6) with 

du /dx = 0 produces Equation (150), with the following values of A and B 
0 

for turbulent flows: 

A = ( 0-8 ) [ 0.028 (1+n) (2+n)- j 0,8 

'n+1 7 n 


(170) 


B = 0.2 


(17D 


The following empirical correlations of the exponent n and factor B have been 
developed to extend the applicability of this result over the range of Reynolds 
numbers and pressure gradient effects discussed by Schlicting (Reference 15, 
Chapter XXI and ] XXII): 


2.183677 + .573308 log R - . 01 36455 (log R) z for R < 1.7x10 


n = 
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7.6 for R > 1 . 7x10 
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1 0.2 - .061 exp (R - 1.7x10 ) 
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fof R < 1.7x10 


for R > 1 .7x10 
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(174) 


where p- = 3p/3x, L Is the reference length used for non-d I mens i ona 1 1 za 1 1 on of 
the inviscid flow, and Ap is the overall static pressure ratio across the stage. 


Blade Wakes 


The main interest in the wakes of a blade row resides in the reaction of the 

blades of a succeeding blade row to traversal of each of the wakes of the first 

row. The reaction is principally that of an airfoil passing through a nonuniform 

20 

inviscid flow field (i.e., a "gust") . However, the degree of nonuniformity 
(and the intensity of reaction thereto) depends on the diffusive properties of 
the viscous flow; the axial distance between blade rows is paramount in this 
regard. Although further diffusion can occur downstream of the plane of the 
leading edges of the second row, the identity of the individual wakes of the 
first blade row will be lost to a large extent after impingement upon the second 

row. Therefore, in the present program the inviscid slipstreams and viscous 

wakes of the first blade rows are only calculated up to the plane of the leading 

edges of the second row, ?.e., the interface of domains k and 5- At this sta- 

tion the inviscid flow field at the entrance to domain 5 at any Instant is con- 
sidered to be the composite of the inviscid and viscous solutions at the exit 
of domain k. In other words, the identity of the individual wakes of the first 
blade row is lost downstream of this station, but their momentum and energy 
defects are transferred to the inviscid flow field through definition of cir- 
cumferential distributions of flow properties at this station which are com- 
posites of the inviscid and viscous wake solutions. The manner in which the 
composite solution is formed is described following the development of the wake 
solut ion. 

6 1 



In view of the intended application of the wake sol ut ion, g reater emphasis 
has been placed on accurate description of the velocity defect within the wakes 
than is given to description of the boundary layer velocity profile. Accordingly, 
a more rigorous theoretical model of the wake has been derived, although it still 
lies within the framework of a quasi-steady integral method of analysis. 


Governing Equations and Boundary Conditions 

The standard boundary layer equations are assumed to apply along the blade 
slipstreams. The x, y coordinate system defined in connection with analysis of 
the slipstream (Figure 7) forms the present boundary layer coordinate system 
which is shown in Figure (11). The governing equations are, therefore, stated 
as : 
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H = h + ~ (u 2 + v 2 - tt 2 r 2 ) (179) 

All time derivatives have been neglected and the standard boundary layer ap- 
proximations have been invoked, (e.g., v << u) . The Prandtl number y C /k 

P 

has been assumed equal to unity. The viscosity coefficient y can be considered 


*ln this context y denotes an outward normal from the slipstream; thus y = |y|. 
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FIGURE II. SCHEMATIC OF WAKE COORDINATE SYSTEM 


to represent either the laminar value associated with molecular diffusion or 
an effective "eddy vi scos i ty" associated with turbulent motion. in view of 
the extreme remoteness of the possibility that a laminar flow could exist In 
the wake of a fan or compressor blade under normal operating conditions in 
either full-scale flight or a sub-scale test facility, only the turbulent 
case is considered herein. In this connection, it should be emphasized that 
the previously noted comparisons of inviscid and viscous time scales should 
be extended to_ include the character i st i c frequency range of turbulent fluc- 
tuations. The present usage of "eddy viscosity" concepts and associated re- 
sults (der i ved, f rom experimental observations in steady inviscid flows) is 
predicated on the assumption that the frequency range of inviscid flow oscil- 
lations induced by blade row interactions does not overlap the turbulence 
frequency range. Consequently, the two unsteady flow phenomena can be un- 
coupled. This assumption is of course suspect, but its removal requires re- 
search into the fundamentals of energy exchange between the inviscid flow 
and turbulent motion that is well beyond the scope of the present study. 

The boundary conditions pertaining to the wakes are summarized as: 


At y = 0: 

v = 
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( 1 80 ) 
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At y y^.: 

u = 

(u 2 + v 2 )~ 

( 182 ) 
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At x = x^. : 
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^boundary layer 

(185) 
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The subscript f denotes the instantaneous inviscid flow values measured at 
some suitably defined front between the inviscid and viscous flows, which im- 
plies y^. < 2-nr/H. The boundary conditions for the case when the wakes merge 
to form a completely viscous flow will be stated in a following paragraph. 


The conditions expressed by Equation (l8l) are approximate, since the 
wakes will not, in general, develop symmetrically about the slipstream posi- 
tion defined by y = 0. However, the displacement of the line of minimum ve- 
locity and zero normal gradient of velocity from the slipstream position will 
be neglected. Furthermore, the asymmetric character of the wake, due to the 
differing boundary conditions which may exist on each side of the wake as 
y ■+ y^., as well as the asymmetry of the boundary layer solutions of the blade 
trailing edge, will be approximated by using two symmetric solutions; one per- 
taining to 0 < y < y^ and the other to 0 < y < y^. 

At sufficient distance downstream of the blade row the wakes will entrain 
the entire inviscid flow field and merge to form a fully turbulent flow. In 
this event, the boundary conditions given by Equations (182), ( 1 8 3 ) and ( 1 84) 
are replaced by: 


At y = y 
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The analyses pertaining to independent wakes and merged wakes differ correspond- 
ingly as outlined below. 


Analysis of Independent Wakes 
2 1 

The method of integral relations is employed using a single strip to cover 
the region from the wake axis (y = 0) to the wake front (outer edge) (y = y^ ^ 
or y f ). This approach is equivalent to use of the Karman momentum integral 
equation for the boundary layer; however, the surface boundary conditions are 
replaced by a set of differential equations corresponding to Equations (176) 
and (177) evaluated at y = 0, where Equations (180) and ( 1 8 1 ) pertain: 
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(The subscript m denotes minimum value, which in accord with the assumption 
of symmetry occurs at y = 0.) 

The integral form of the continuity equation results from multiplication 
of Equation (175) by dy and integration from y = 0 to y = y^(x) (where y^. 
denotes both y^ ^ and y^ u ) Observing the above noted boundary conditions at 
these limits, integration yields: 

*f _ *f 
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Corresponding integration of Equations (176) and (177) and substitution of 
Equation (191) into the resulting equation yields; 
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(192) 
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Equations ( 1 9 1 ) , (192) and (193) each represent a pair of equations; one for 

y f = y f,r u f,n- v f = v f,r etc - and the other for y f = y f,u* u f = u f,u' etc - 


The y coordinate is transformed according to: 
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and transformed wake thicknesses are defined by: 
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The eddy viscosity coefficient is assumed to transform according to : 


f ,£ 


f ,u 


(196) 


where e. denotes an equivalent incompressible value. The above transforma- 
tions reduce the system consisting of Equations ( 1 89 ) through (193) to that 
of an equivalent incompressible flow. The eddy viscosity in incompressible 

JL. 

flows has been experimentally correlated by : 


*For the sake of clarity, discussion of the basis of the selected eddy viscosity 
law and velocity profile representations is deferred to the end of this derivation. 
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e. = .032 (u f - u ) 6 

i t m 

In view of the noted asymmetry in the present application, Equation (l37a) 
is interpreted as: 


e. = .032 ((u, „ + u, )/2 - u ) (6„ + 6 ) /2 

1 f , £ f ,u m £ u 


The velocity profile is represented by a corresponding approximation: 

/ _ 


“f ‘ = { 
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1 (u £ - u ) (1 + cosirn / 6 )/2 for 0 < y < y £ 
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A similar representation is assumed for the relative total enthalpy profile, 
however, the possible existence of a different thickness for the relative 
total enthalpy defect than for the velocity defect is recognized: 


H f - H = 


f (H f,£ • V (1 + cos, V 6 T,i )/2 for 0 < y < y T> 


(H - H ) (1 :+ cos Tin /6-r )/2 for 0 < y < y_ 

f,u m u T,u — — T,u 


where 6_ f <5 „ and 6 T ^ 6 is assumed. In this connection, it should be 
T, £ £ T,u u 

pointed out that compressor or fan blades are usually uncooled and, therefore, 

operate at a wall temperature corresponding to the average absolute total 

temperature (or total enthalpy) of the airstream in which they are immersed. 

Thus, the difference H_ - H will be on the order of the difference between 

t m 

the local, i ns tan taneous total enthalpy and the average value experienced by 
the blade. !t is anticipated that in general: 


H f - H « H f 


Furthermore, it is well established that 
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(200) 


68 


( 201 ) 


u << u, 


prevails within a turbulent wake a very short distance from the trailing 
edge. 1 1* therefore, foi lows that: 


p/p. 


1 


Also in connection with Equations (198) and (199), it should be noted 
that double-valued second derivatives of velocity and total enthalpy at 
y = 0 result from the two forms of each equation. Therefore, the second 
derivatives at y = 0 are approximated by: 


a 2 " 

(^) 

3n 


2 TT 


(5 j + V 


■n ((u f + u f )/2 - u ) 

2 f,u m 


(3V) 

9n 2 


2tt 


^T,£ + 5 T,u^ 


■ 9 (( H, + H, )/2 - H ) 

2 f,£ f,u m 


Substitution of the above-stated profile representations, eddy-viscosity 
law, and transformations (e.g.. Equation (178)) into the system given by Equa- 
tions (189), (190), (192) and (193), and using p/p^ * 1 In evaluation of the 
Integrals, yields the following set of ordinary differential equations: 


du 


dx 


2 

1 dp .032 tt r /- , — \ /n - , 

- — + 77 7 j - r [( U 4T „ + tl r ,,)/2 “ u 1 

P (6 . + 6 J f , 2. f ,U m 

m dx i u 


. r , 2 dr 

+ ft r -7— cose 
dx 


dH^ _ -032, («,+«„) 

d * (6 T.!. + 6 T,u )2 


+ “f,u )/2 ' “m 1 [(H f,* 


+ H, )/2 - H ] 
f,u m 


(202) 


(203) 


(20k) 


(205) 


( 206 ) 
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, d6 ; 


1 


du du* . 

( — - — -^) - (- 


*7 dx (u^ . - u m ) dx dx 


2 du f,i , 1 dp f,i 


u, . dx 

f , ' 


f , i dx 


cos dA _ 2Q r cosH; d r ? for I=£ and u 


A dx 


U,- . (U, . “ U ) 

f , i f , t m 


dx 


1 d6 T * 

1 T, t 


dH dH, . 

r m f , I 


du r . 


T , i dx 


(hU . - H ) dx dx 
f , i m 


) - (: 


dp. 




u £ . dx 
f , i 


f , i dx 


, dH r . 5. (u. . - u ) 

1 / f,i^ \ f,i m cos d> dA * _ . A 

— ( — ) -r ^ ^ - — t— — , for i = Z and u 

- o ■ i \ A dx* 


u r . dx 

f , i 


T . (H, . - H ) 
T , i f , i m 


A more specific definition of the frontal locations ^ and u is now 
required to carry out the numerical integration of the above system. The 
frontal location can be unambiguously defined by matching the mass flow con- 
tained in the wake at any streamwise position to that in a corresponding in- 
viscid streamline, namely: 


m = 


(p u dy) . . . , 

i nv i sc i d 


y . 

[ 


(p “ dy) wake “ °f "f y f 


Thus ; 


whe re : 


P f u f 


(p u dy) . . . , 

inviscid 


u r = u. 
f i n v i sc 


id ( *- y = v 4) 


P • • • J (x, y = Yr * t) 

inviscid f 


etc . 


( 207 ) 


(208) 


( 208 ) 

( 210 ) 

(211a) 

(211b) 
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and : 


du, 


dx 


3u. . . ,1 

i nv i sc i d j ( 

p u invi sc i d 

dy f 

(212) 

. 3x J „ 

[ay 

| - - - ^ dx 


x,y=y f ,t 

1 x,y=y f ,t 



The x derivatives of and must be correspondingly defined. Since p/p^ « 1; 


dy £ . dS. 
r f , i _ i 


for i = £ and u 

dx dx 

Therefore, Equation (207) must be rewritten as: 


(213) 


1 


dfi. 


dx 


. - u ) 
f , i m 


du 


dx 


<%) 

3x 


u r . 3x 

f * • 


3u _J 3pv _ cos <j> dA _ 2fT rcostj) 

P r t "i T. A dx 


dr 


f , i 3x 


u r . (u r . - u ) dx 
f , i f , i m 


. 3u^. . - 2u 

J_ . f > > m 

< 5 . 


(liL) + _L 12. 


u f i ( u f ; ■ u m ) ay 

T , l I y i m 


f 3y 


( 214 ) 


for i = £ and u 

where 3/3x and 3/3y denote the inviscid flow values at (x, y = y^., t) as used 
? n Equat ion (212) . 


To within the approximation contained in the above system of equations, the 
displacement, momentum defect and thermal energy defect thicknesses in the wake 
are defined by: 


u, . - u 
f , i m 



1 

i 

2u . 








u . - 


0 _ 

*0 

11 

f , i 




2U r . 

f , 1 


0 

<o 

1! 

h' . 
f , 1 


T,i T , i 


2h m 


for i = £ and u 


for i = £ and u 


for i = £ and u 


(215) 


( 216 ) 


(217) 
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Thus, 6j = 0. in this context. However, in the boundary layer at the trailing 
edge 6* ^ 0; therefore, the wake solution cannot match both the displacement 
and momentum thicknesses of the boundary layer solution, and a mismatch in one 
or the other must be accepted. Since the displacement thickness is to be used 

JL 

as an effective boundary of the inviscid flow, a mismatch in 6 is considered 
unacceptable. Therefore, the initial conditions for the wake solution are 
taken to be: 

u (0) = 0.0 

m 

h' (0) = (C T) - ~ Si 2 r 2 

m P blade 2 

trailing edge 

5 (0) = 2 6 W , . , , 

l blade upper surface 

.trailing edge 

6 ( 0 ) = 2 6 '' U1 . . * 

u blade lower surface 

trail in g edge 


11 

o 

1- 

«* (o) 

(221a) 

T,u<°>= 

5 u (0) 

(221b) 


(The initial thermal defect thickness is equated to the initial velocity de- 
fect thickness ?n lieu of any more precise information. Should a more de- 
tailed boundary layer solution be incorporated in the model at a later date, 
the initial wake conditions can be modified accordingly.) 

Analysis of Merged V/akes 

In the foregoing discussion, it has been implied that the wake of each 
blade grows into a known inviscid flow field. However, at a sufficient dis- 
tance downstream of the blade row the wakes will entrain the entire inviscid 
field and merge to form a completely turbulent field. In this case the pre- 
viously stated analysis must be modified accordingly. Following the previously 



(218) 

(219) 

(220a) 

(220b) 
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outlined procedures, using the same small perturbation approximations, and 
adopting the same profile representations , the integrated form of the con- 
tinuity equat ion becomes: . 


- — (u + u f ) = 
dx . m . f 


Cu + u,)' SZL± dA 

m t A dx 


where now u f = = u f and 6 = (6 + 6 )/2. Evaluating the momentum eq.ua- 

lljA/ljU XU _ 


tion on the axis (y = 0) and at y = y f £ = 2-iTr/N - y f u gi 


ves v 


du 


dx 


J_ d£ + 


2 2 
7 r 


m dx 


— 0 r\r* ' 

032 (u r - u ) — + fi r -j— cos<J> 

i rn ■ « » d x 

Zo 


and 


du , 


— — - -032 (u f - u ) — + fi 2 r 4^- cos<{> 

f m 25 dx 


dx f dx 

Substituting from Equation (22k) into (222) gives: 


du 


dx 


h ¥ + - 032 < G f - =/ ^ 


f dx 


26 


r , - \ cos <j> dA 
u f u m u f dx 


_ 2 dr 
fi r -j— cosi 
dx 


Adding Equations (223) and (225), and noting that - p m << p^» yields: 

du 


(u f - u ) 
t m 


dx 


2 2 

_y>/ - - \ tt - t - \ cos 4> dA 

.064 (u r - u ) — - u, u + u.) — t — *- -t~ 

f m f m r A dx 

2o 


Thus, the rate of chanqes of u c and of u are determined from Equations (222) 
5 3 r m 

and (226). Furthermore, the pressure is determined from the sum of Equations 
(223) and (224) : 


. . . du 

(_L + J_) dp = - (u 

p f p m dx m dx 


m - du f\ 0 dr 

+ u„ -) + 2fi r ^ cos(J) 


f d ; 


-Note that the cosine function is particularly appropriate for this problem 
due to the periodic character of the wake. See Reference (15), Chapter 
23, p. 604. 


(222) 


(223) 


(224) 


(225) 


( 226 ) 


(227) 
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Evaluation of the energy equation on y = 0 and on y = y^ = 2irr/N - y^ u 
gives: 


dH _ 2 

i _JH = .032 (u f - u ) (H £ - H ) — 

m dJ f m ' m 26 


( 228 ) 


and 


dH , 


dx 


.032 (u f * u) (H f - H ) — 

f m f m 26 


(229) 


from which the rates of change of H and H r can be obtained. Note that 

m 1 

5 = 6 is required downstream of the point of merger. 


it is remarked parenthetically that in the speciai case of o’A/dx = 0, 

Equations (222) and (226) can be combined to yield the well known result 

that the maximum velocity differential (u f - u ) must asymptotically decay 

t m 

linearly with increasing distance, in contrast to asymptotic square root 

15 

decay rate for the individual blade wakes. 


Eddy Viscosity Law and Velocity 
Profile Representations 

A variety of formulas have been proposed in various investigations to re- 
late the eddv viscosity defined bv t h e Reynolds stress, namely: 


u v _ . 

(3u/3y) k t — fcclen: 


(230) 


to the mean properties of a turbulent flew. Tbe present formulation will draw 
upon one of the earliest and most generally accepted eddy viscosity laws; 
however, its validity can only be established by comparison with experiment and 
it should always be regarded as tentative and subject to revision. 


A "universal" number characterizing the eddy viscosity in a turbulent, in- 

23 . 

compressible wake was derived bv Townsend " based on equilibrium of energy 
among the large-scale eddies: 



constant 


(231) 


R_ = (u r - u ) £ /e - 
T f m o 


where £ Isa wake width parameter defined by the requirement that the ve- 
o 

Tocity profile conform with: 

u, - u = (u_ - u ) exp (-y 2 /2£ 2 ) 

f f m o 

2 3 

Townsend obtained the value R^ - 12.5 from a correlation of his data on the 
wake of a cylindrical rod in the form: 


_ -2 

- u = (u f - u ) exp {- ~ — i 

f m 2 C d (x - x )d 


from which £ can be identified as: 
o 


£ = (C , (x - x )d//2rr" R ) 

o d o T 


where C, is the drag coefficient, d the rod diameter and x an effective 
d 3 o 


origin 


A somewhat more useful definition of £ can be obtained as follows. 

o 


We note that Townsend's wake data can also be accurately correlated in terms 

2k 

of Coles wake function as indicated in Figure (12): 


u - u 


u r - u 
f m 


W(n) 


where n = y/6 and 6 is the value of y where u ~ u^. W(n) is a tabulated func* 
tion based on correlation of turbulent boundary layer data, and is known to 
have wide applicability to compressible as well as Incompressible flow through 


suitable transformation 


25 


W(n) can be accurately approximated (see Figure 12) 


by W = 1 - cosTrr). Therefore, 


= ( u f " u m ) 
t m 


(1 + costth) 


The above velocity distribution is matched to the desired exponential form. 
Equation (232), by requiring that: 
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( 232 ) 


(233) 


(23*0 


(235) 


( 236 ) 





I 


exp {- y 2 / 22 , 2 ) dy = y 


(J + cos-rry/S) dy 


2, = 6 /'/ 2 t \ ^ O.AO 6 

o 

Therefore, if the velocity profile is represented by Equation (23&) , the 
eddy viscosity is consistently represented/ by : ... 

e = 0.032 (u f - u m ) 6 

Parenthetically, it is noted that in terms of the exponential form, the above 

value of 2, gives u = .957 u, + .0^3 u at y = 5. 
o t m 

it is emphasized that the above results pertain to incompressible flows. 

The extension to compressible flows by transformation has been discussed by 
22 

Ting and Libby . In particular, the postulate that the shear stress remains 
invariant under density transformation (see Reference 22) has led to satis- 
factory results at supersonic and hypersonic conditions and, thus, should be 
quite accurate for the presently considered transonic and low supersonic flows. 


Composite Solution 

As previously indicated, the viscous and inviscid solutions obtained at the 
exit of domain A are used to define a compos i te i nv i sc i d solution for the en- 
trance station of domain 5- The composite solution is defined as follows: 


(ucosb) . 

viscous 


v t scous 


inviscid 


i nvi sc i d 


a < 0 


0 < a < 1 


a > 1 


+ V 


a < 0 


v = < 


(us i n<j>) 


viscous inviscid 


(l-a) (usini) . + v. . . , 

viscous inviscid 


v . 


i nvi sci d 


H . 

VI scous 


= s ( 1 -a) H . + a H . . . , 

] vi scous inviscid 


H. 


Inviscid 


P = P • • • j 

r inviscid 

where : 

a = (y -6*)/(<5-<T) 
y = I y-y - I COS^ i = £ and u 


0 < a < 1 r 
a > 1 


( 241 ) 


a 0 i 

0 < a < 1 l 

a > 1 I 


(242) 


(243) 


(244) 

(245) 


The variables denoted by u . and H . are obtained from the pro- 

v i scous v i scous 

file representation given by Equations (198) and (199) using the viscous solu- 
tion u , H , 6. and 6_ . and the frontal values u r . and H r . ( i = £ and u) . 

The range of values of y is selected to match the undisplaced y coordinates 
of the grid points on those two grid columns of domain 4 which bound the first 
grid column of domain 5, and composite solutions are generated at these two sta- 
tions, Equations (240) - (243). The inviscid solution for the first grid 
column of domain 5 is then obtained by linear interpolation of the composite 
solutions, for use in accord with the procedures outlined in the sub-section 
entitled "Periodicity Condition". Thus, the composite solutions simply replace 
the inviscid solution which would be used to carry out the interfacing of domains 
4 and 5 in the absence of viscous effects. 


*The transformed ci rcumferent ia 1 coordinate v used for the inviscid solution in 
domain 4 (and also domains 3, 5, 8 and 7) obviously includes the local displace- 
ment thickness in definition of the computational boundaries; I . e . , in Equation 
( 37) y j = r0 j +5 •" (i = & , u) . Thus, the displaced computational domain does not 
span the entire blade-to-blade passage, whereas the undisplaced domain does. 
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RESULTS 


The method of analysis described above has been implemented in a FORTRAN 
26 

code and tested with respect to several single stage transonic fans for which 

27 28 

experimental data is available ’ . A combination of design information and 

measured data for these fans has been used to define the initial and boundary 
conditions. The "infinite duct" inlet and discharge boundary conditions were 
used in the calculations reported here; limited tests with the "open end" con- 
ditions were also conducted, but detailed comparisons with the "infinite duct" 

model were not carried out. The meridional plane analysis of Katsanis and 
29 

McNally was used to define the streamsheet thickness and radius by tracing a 
selected streamtube through the stage. 

As a prelude to analysis of rotor-stator interactions, several isolated 

rotor analyses and an isolated stator analysis were carried out. The stator 

analysis will be discussed first, since comparison of the present method and 

30 

the relaxation method of Katsanis is possible for this case, which offers a 

direct assessment of the program accuracy, independent of the viscous and three- 

dimensional effects present in the transonic rotor data. Since the stator is 

subcritical at the selected operating condition, the velocity-gradient tech- 
30 

ntque provides a reasonably accurate benchwork solution for a shock-free 

27 28 

flow. The experimental data * for the rotors will be used primarily to 
indicate the predictive capability of the present method with respect to the 
structure of passage shock systems. 

o *7 

The stator is part of a 1 500 fps (^57 m/s) tip speed transonic fan stage' , 
which will be described more completely later in connection with the rotor-stator 
interaction analysis. The stator row has bG blades and a hub to tip ratio of 
0.6. The selected operating point corresponds to Reading 137, of Reference (27) 
which is an open throttle, 100% speed condition. The selected streamsheet fel- 
lows the casing wall and has a thickness of 1/3 of \% of the casing radius at 
the inlet station (one chord length upstream of the rotor). The inlet Mach num- 
ber is 0.6 and the inlet flow angle is -30.7°. The stator is intended to produce 
a purely axial exit flow. The blade section at the casing is shown at the top 
of Figure (11). A grid network consisting of 17 columns (axially) and 9 rows 
(circumferentially) in each of three domains (i.e., ^59 points, excluding external 
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grid points) was used for this case, as well as the isolated rotor cases to be 
discussed next. The solutions converged within 10^ time steps with this grid 
point density, and each required approximately three minutes execution time on 
a CDC 7600 computer. 

An inlet Mach number of O .65 and an inlet flow angle of -30.0°, with neg- 
ligible circumferential variations, were obtained from the present stator so- 
lution. The outlet flow angle varied circumferentially from 6.6 to 7-0°. The 
relaxation solution was carried out for the same inlet Mach number of 0.65, but 
using the design inlet and exit flow angles of -30.7° and 0.0°, respectively. 

(Specification of the exit flow angle replaces the trailing edge Kutta condition 

30 

in the relaxation method .) Sixty-one grid columns and 20 grid rows were used 
in the relaxation solution. The surface Mach number distributions are shown on 
the bottom of Figure (13)- It can be seen that the effect of nose bluntness 
(which is included in the relaxatign solution^ 0 ) is essentially confined to t he 
first 5 - 10% of the chord. The compression surface solutions are in quite good 
agreement up to the last 20% of the chord, where the effect of the manner of en- 
forcing the Kutta condition is evident. The suction surface solutions exhibit a 
fairly uniform difference; the present method results in a somewhat higher Mach 
number over almost all the surface. In view of the doubled grid point density 
in the ci rcumferential direction used in the relaxation solution, it must be re- 
garded as numerically more accurate; however, the accuracy of the velocity- 
gradient approximation for a local Mach number approaching unity is uncertain. 

27 

The rotor for this stage has 44 blades and a hub-to-tip ratio of 0.5. The 
tip speed is 1500 f ps , (457 m/s), producing an inlet relative Mach number of 1.526 
at the design point. The tip diameter is 36.5 inches (0.927 m) and has an axial 
chord length of 1.7 inches (0.0432 m) at the tip. The rotor has a shroud (vibra- 
tion damper) located about 40% of span in from the tip. The stage also included 
2k variable camber inlet guide vanes, located slightly more than one rotor chord 
length upstream of the rotor. Under the presently considered conditions the 
guide vanes were set to zero camber, to produce a purely axial inlet flow to the 
rotor. The same operating point (Reading 137) was examined, at which the stage 
total pressure ratio was 1.48 arid the inlet relative Mach number at the tip was 
1.49. The streamsheet was again assumed to be a very narrow layer along the 
casing wall (having a thickness about 1/3 of 1% of the radius). 
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SURFACE COORDINATE, rfl/c 




FIGURE 13. STATOR BLADE TIP SECTION GEOMETRY AND COMPARISON 
OF SURFACE MACH NUMBER DISTRIBUTIONS 
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An inlet Mach number of 1.47 was obtained in this case, and the rotor total 
pressure ratio varied (circumferentially) from 1.37 to 1.50, in good agreement 
with the design conditions. The rotor pressure contours obtained from the pre- 
sent solution are displayed in Figure (l4a) and compare favorably with those ob- 
tained from arrays of fast-response pressure gages in the casing wall , shown 
in Figure (14b). A more direct comparison is offered in Figure ( 1 5 ) which pre- 
sents the blade surface pressure distributions. The differences near the lead- 
ing edge are probably attributable to the effects of nose bluntness. The ob- 
served regions of compression-expansion preceding the passage shock (z/c = 0.7 
on the compression surface and z/c - 1.0 on the suction surface) are indicative 
of local boundary layer separation bubbles. The part-span shroud may also be 
generating a shock wave which interacts with the blade shock system. In view 
of these factors, the general agreement between theory and experiment is re- 
garded as quite satisfactory. 

The three-dimensional shock structure near the tip of a similar transonic 
rotor was visualized using pulsed laser holography in Reference (28). In this 
example, the rotor has a tip speed of 1600 fps, (488 m/s), and a design pressure 
ratio (for the stage) of about 1.5. The tip diameter of the fan is 28.74 inches 
(0.730 m) and the (axial) chord length is 1.804 inches (.0458 m) . The inlet hub- 
to-tip radius ratio is 0.46 and the nominal inlet relative Mach number is 1.6 at 
the tip. The rotor has 40 blades, and includes a midspan vibration damper lo- 
cated about 30 % of the span in from the tip. 

The streamsheet thickness and radius distributions were obtained by traci.no 
a streamtube located about 5% of the span in from the tip and having 0.11 inches 
(.00279 m) thickness at the inlet station one chord upstream of the blade row. 

Since its thickness is less than 1% of its radius, the streamsheet radius closely 
follows the casing of the machine. Three operating conditions were examined; 

(a) the Resign point (Reading 128), (b) 100% design speed with a pressure ratio 

of 1.7 (Reading 126) and (c) 90% speed with a pressure ratio of 1.5 (Reading 106). 

The computed pressure contour map for the design case is shown in Figure ( 16 ) . 
The presence of a weak shock off the blade leading edge is evident; it subsequently 
reflects o p f the lower blade at about the 85% chord position and then apparently 
merges into a stronger shock which crosses the Dsssage from the trailing edge of 
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CIRCUMFERENTIAL DISTANCE, r9 (IN.) 



FIGURE 14a CALCULATED STATIC PRESSURE CONTOURS AT ROTOR TIP 
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CONTOUR PRESSURE 
NUMBER (PSIA) (P/Pcq) 

1 9 .7647 

2 10 .8496 

3 1 1 .9346 

4 12 1.0195 

5 13 1.1045 

6 14 1.1895 

7 15 1.2744 

8 16 1.3594 

9 |7 1.4444 

10 18 1.5293 

11 19 1.6143 

12 20 1,6992 

13 21 1.7842 


AXIAL DISTANCE, Z (IN.) 


FIGURE 14b. MEASURED STATIC PRESSURE CONTOURS AT ROTOR TIP 
(READING 137 FROM REFERENCE 27) 
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FIGURE 15. PRESSURE DISTRIBUTIONS ON ROTOR BLADE TIP 
SURFACES (READING 137). 
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FIGURE 16. COMPUTED ISOBARS FOR READING 128 OF REFERENCE (28) 
M^ = 1.6 AND PRESSURE RATIO = 1.5 



n 


the lower blade to about the 85% chord position on the upper blade. These fea- 

• 28 

tures are generally in accord with the reconstructed holographic observations ; 
however, the interaction between the tip leakage vortex, which stands on the suc- 
tion {i.e., low pressure) surface of the lower blade, apparently causes the rear- 
ward portion of the leading edge shock to bend forward somewhat such, that it re- 
flects off the lower blade at about 75% chord (rather than 85%), and the trailing 

edge shock is correspondingly displaced upstream. The shock off the vibration 

28 

damper also appears to be interacting with the blade shocks . In addition, the 

assumption of a sharp leading edge may have reduced the strength of the leading 

edge shock. The computed total pressure ratio at the discharge station varied 

(circumferentially) from 1.^5 to 1.51, in reasonable agreement with the overall 

2 8 

stage pressure ratio of 1.505 reported at this operating condition . The rela- 
tive Mach number at the inlet converged to 1.51 as compared to the initial value 
of 1.53. 

In the second case (Reading 126), the higher back pressure produced a modest 
forward shift of the trailing edge shock near the upper blade, as can be seen in 
Figure (17). The solution along the suction surface of the lower blade, as well 
as the forward 50-60% of the compression surface of the upper blade, is virtually 
unchanged, since forward propagation of the higher back pressure is terminated at 
the trailing edge shock. No experimental data is available at this operating 
cond i t i on . 

As can be seen from Figure (l8), the reduced wheel speed of the third case 

(Reading 106) results in a lower relative inlet Mach number, i.e., about 1 . ^ , 

and a further shift of the trailing edge shock, to form a normal shock across the 

28 

passage. However, the holographic reconstruction indicates that the rotor is 
"unstarted" at this condition, i.e., the normal shock stands across the passage 
at the leading edge of the upper blade. The solution for this case was perturbed 
several times, e.g., by increasing the back pressure, by altering the streamsheet 
thickness distribution, and by reinitializing the data to an "unstarted" condition, 
in an unsuccessful attempt to produce an "unstarted 11 solution. Therefore, this 
disparity between theory and experiment must be attributed to one or more signifi- 
cant features of the acutal flow field identified above but not incorporated in 
the present analysis: finite nose bluntness, viscous effects, and three- 

dimensionality. Nose bluntness effects are believed to be significant in several 
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FIGURE IE COMPUTED ISOBARS FOR READING 106 OF REFERENCE (28) 
M^=l.4 AND PRESSURE RATIO = 1.5 
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respects, and should be included to enhance the accuracy and expand the range of. 
applicability of the present analysis, but in the present case do not appear to 
be sufficient to cause the passage to "unstart". The same comments can be made 
with respect to viscous effects on the blade surfaces. However, the tip leak- 
age vortex and the vibration damper shock present major three-dimensional per- 
turbations. to the flow near the tip section, which are not easily accommodated 
in a two-dimensional model, but could cause the noted discrepancy between the 
present results and the observed unstarted shock system structure. 

Finally, the rotor-stator interaction has been examined at the open throt- 
tle, 100% speed operating point (Reading 137 of Reference 27) previously selected 
As pointed out earlier, the stage includes 44 rotor blades and 46 stator blades. 
Since the flow leaving the rotor is supersonic (M - 1.1) upstream propagation of 
disturbances from the stator should be cut-off beyond the steady wave front in- 
tersecting the trailing edge of the lower rotor blade of the passage. If the 

flow were uniform in the passage, this wave front would intersect the upper rotor 
blade at about the 40% chord position; due to the nonuniformity of the flow it 
actually intersects the upper blade at between 60 and 80% of chord. (See Figures 

14a and 14b.) The rotor flow field upstream of this wave front converges rapidly 

to an essentially steady solution, whereas the stator is subjected to sequence of 
perturbations which travel both upstream and downstream and, therefore, converge 
to a periodic solution somewhat more slowly. The calculation was carried out for 
one complete revolution of the rotor, which required about 30 minutes of CDC 7600 
computer time, using a grid network consisting of 12 grid columns and 9 grid rows 
in each of 5 domains. Approximately 1/4 of a revolution (500 time steps) was re- 
quired to achieve an asymptotic solution in the rotor passage, and about 1/2 rev- 
olution (1000 time steps) should have been sufficient to attain an asymptotic 
solution in the stator. Unfortunately, a minor coding error in the application 
of the phased boundary conditions was not discovered until 3/4 of a revolution 
(1600 time steps) had been computed. Upon correction, a periodic solution was 
achieved during the last 1/4 of the revolution. The entire solution was not 
repeated due to the computing time requirement. 

Provision to integrate the surface pressure distributions to obtain normal 
force and moment coefficients was not included in the code, making display of 
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the periodic solution rather cumbersome. Therefore, the pressure difference 
across the blade at the midchord point has been selected to represent the blade 
force. The midchord pressure differential across the stator is shown in Figure 
(19) for a period covering approximately the last 1/A of the rotor revolution. 
Presence of a periodic solution formed by superposition of two very distinct 
waves is apparent. One wave appears to have a peak-to-peak frequency twice the 
rotor passing frequency, and the other has a peak-to-peak frequency twice the 
stator passing frequency (relative to the rotor). Fourier analysis of the solu- 
tion over the period of 1/2 a revolution would be required to quantify the spec- 
tral components of the blade row interaction field. 

As indicated above, the asymptotic flow field upstream of about the mid- 
chord position of the rotor should be steady in the rotating frame of reference 
of the rotor. However, as can be seen in Figure (2o)» a periodic pressure dif- 
ferential across the rotor midchord location is obtained, and in fact even the 
Inlet pressure distribution exhibits periodicity in the rotating coordinate sys- 
tem. However, the amp] i tide of these oscillations is much smaller than those 
found in the stator. in view of the noted supersonic character of the relative 
fl ow through the rotor, these oscillations must be attributed to numerical propa- 
gation. Once the oscillations produced by the stator reach the leading edge of 
the rotor they are correctly able to propagate upstream, since the axial velocity 
component is subsonic. This spurious numerical propagation could be eliminated 
by observing the correct domain of dependence in formation cf the difference 

operator, along the lines of the type-dependent difference operator used in the 

2 

relaxation method developed by Murman and Cole . 

SUMMARY AND RECOMMENDAT i ONS 

A numerical method of solution for inviscid, transonic flow through a set of 
interacting cascades has been described in detail. Particular attention has been 
devoted to the statement and method of implementation of boundary conditions. The 
solution algorithms employed at the interior and boundary points of the computa- 
tional grid have been described in detail. 
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FIGURE 19 STATOR MID-CHORD PRESSURE DIFFERENTIAL DURING 
QUARTER REVOLUTION OF ROTOR 
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FIGURE 20- PRESSURE VARIATIONS AT ROTOR MID-CHORD AND AT 
INLET DURING QUARTER REVOLUTION 
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Several numerical examples have been carried out to demonstrate the cap- 
abilities and limitations of the present method. For example, a sharp blade 
leading edge is assumed in the present mode 1 . Com pa rison with a rel axat i on 
method solution for a stator blade row having blades with small but finite 
leading edge radius and operating under subcritical, transonic flow conditions 
indicates that the effect of nose bluntness is felt over the first 5 to 10 % of 
the chord, which represents 10 to 20 nose radii for this case. In addition, 
comparison with experimental transonic rotor data indicates that the leading 
edge shock is generally stronger than that obtained from the present analysis, 
which is attributable to neglect of the nose radius. Thus, the range of ap- 
plicability and degree of accuracy of the present computational model could be 
enhanced by inclusion of finite nose radius. However, the accuracy of the solu- 
tion over the major portion of the chord length, under subcritical flow condi- 
tions, was found to be quite good with only 9 grid rows used from b 1 ade- to-b 1 ade ; 
obviously it could be further improved by a finer grid network. 

Examination of several transonic rotor solutions for "started" operating 
conditions indicates generally good agreement with the overall shock system struc- 
ture; however, boundary layer separation effects on the blade surface, evident 
in the data, produce local departures from the predicted shock structure. A sig- 
nificant disparity between theory and experiment was encountered at a high back 
pressure operating condition. In this case, the present analysis produced a 
normal shock across the passage near the trailing edge of the lower blade, where- 
as holographic visualization showed a normal shock at the leading edge of the 
upper blade, i.e., the rotor was in an "unstarted" operation mode. This condi- 
tion has been attributed to one or more observed effects of the three- 
dimensionality of the actual flow field, namely: presence of a shock off the 

midspan vibration damper, and of the blade tip vortex. Although not explicitly 
observed, separation of the end wall boundary layer may also contribute to the 
disagreement between the present two-dimensional analysis and the experimental 
observations in this case. 

A transonic rotoi — stator interaction case was carried out for one full rev- 
olution of the rotor. Unfortunately, any conclusions regarding the rate of con- 
vergence to a periodic solution were compromised by the presence of a minor cod- 
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ing error which was not detected during the first 3/4 of the rotor revolution. 

A periodic solution was obtained in the stator during the final 1/4 revolution, 
upon correction of the error. The rotor solution became asymptotic within the 
first 1/4 revolution, in accord with the rate of convergence of the isolated 
rotor cases. Integration of the instantaneous surface pressure distributions 
to obtain normal force and moment coefficients would facilitate interpretation 
of the results and is, therefore, recommended as well as Fourier analysis of 
these coefficients to identify attainment of an asymptotic state and to chai — 
acterize the spectral content of their temporal variation. Additionally, nu- 
merical propagation of periodic disturbances through a supersonic portion of 
the rotor passage, which on the basis of the mathematical zones of influence 
should have been inaccessible to upstream travelling waves, has been noted. Al- 
teration of the finite difference operators to eliminate or at least minimize 
spurious upstream propagation through supersonic zones by correctly observing 
the mathematical domains of dependence of the grid points is also suggested. 

Testing of the influence of the blade boundary layer and wake models and 
of the acoustic far field boundary conditions was not included in the present 
investigation due to time and budgetary limitations. In view of the above men- 
tioned areas in which further improvement of the inviscid analysis is considered 
warranted, continued developmental work appears necessary to fully achieve the 
predictive capability of the present method of analysis. Determination of vis- 
cous and acoustic far field effects may logically be regarded as part of this 
next level of development. 
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APPENDIX 


TREATMENT OF VIRTUAL BOUNDARIES . 

i . ■ 

i Periodic Boundary Conditions 

In order to clearly describe application of the periodic boundary con- 
ditions in the present analysis it is useful to illustrate the discussion by 
introducing several rotating devices that demonstrate certain general con- 
cepts with minimal complexity of the geometric configuration. First, con- 
sider a symmetric three bladed rotor developed into a cascade as shown in 
Fi gure (A1 ) : 

CASCADE VIEW FRONT VIEW 



FIGURE Al. ISOLATED BLADE ROW 

Let virtual boundaries of the b 1 ade- to-b 1 ade passages be drawn from the blade 
leading and trailing edges to ±°°. In a reference frame fixed to the blades it 

*These sufaces are not, in general, coincident with streamlines; flow can cross 
them and disturbances can propagate along and across them. 
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is clear that the flow field in each of the three passages must be identical 
if the blades are identical and the inlet and discharge boundary conditions 
are uniform and steady. Therefore, at any instant the flow is circumferen- 
tially periodic. The flow along line 6 is identical to the flow at the 
spatially equivalent points on line ft. The flow along lines a and y is also 
obviously equivalent at any i nstant# The same. argument also applies to points 
in the flow field downstream of the blades. 

Next ; cons'i der the system composed of a rotor and stator each having the 
same number of blades; specifically, consider 3 blades each as shown in Figure 
(A2). 

When the blades are aligned, say at time t Q ; as shown in Figure (A2a) , 
it is clear that, because symmetry produces identical flow channels, the 
flow along and through the virtual boundaries extending from x and between 
leading and trailing edges is the same in each "passage"; i .e. , all influences 
with respect to corresponding points in the respective "passages" are neces- 
sarily the same. For a different relative position, say at some later time 
t + At for example, see Figure (A2b) , the noted symmetry persists since the 
geometry of respective "passages" is the same and again the flow through the 
respective control volumes is identical. The same conclusion as above is 
also reached regarding the equivalence of lines a with y and 6 with ft. The 
same argument again can be applied to the flow field between the two blade 
rows and also downstream of the second blade row. 

It is noted that due to the repetition of identical configurations with 
identical flows it is clear that only one passage of each blade row need be 
computed in order to determine the flow field in the whole periphery of the 
stage. As a consequence of the equivalence of the corresponding interior and 
exterior grid rows, it is also clear that the boundary values, say on exterior 
line <5, can be specified by equating them' to the currently computed values at 
corresponding interior points, in this case on line ft, to enforce the periodic 
boundary conditions for the passage being computed. The same reasoning is ap- 
plied to lines y and a, and. to corresponding virtual boundaries between and 
downstream of the blades. 
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FIGURE A2, ILLUSTRATION OF CYCLIC ALGORITHM FOR STAGE WITH EQUAL NUMBER OF BLADES 
IN STATOR AND ROTOR (N. - 3. N-- 3 ) 



It is not obvious that the case of an unequal number of blades in the 
stator and rotor can be treated analogously; however, the following detailed 
explanation is intended to demonstrate clearly that the above techniques for 
specifying boundary conditions on a single passage can be extended to this 
case by introduction of a phase shift. In this regard a parr of blade rows, 
shown in Figure (A3a) and (A3b) portray the conf i gurat ion of a 3 bladed 
stator followed by a 4 bladed rotor which wi 11 be used to demonstrate the 
unequal spacing techniques. 

The underlying assumption here is that the inlet and discharge boundary 
conditions applied upstream and downstream of the stage are uniform and steady, 
and that the non-steady flow in the vicinity of the stage is stably periodic 
in time. The first question that arises is whether the starting process as- 
sociated with an arbitrary set of initial conditions can lead to an asymptot- 
ically periodic solution through use of the above-described numerical procedure 
for treating spatially periodic boundary conditions. It is useful in this re- 
gard to pictorially describe the physical starting process for a system of blade 
rows set into motion impulsively. This is done for the conf iguration introduced 
above with both the rotational speed and free stream velocity subsonic. 

Figures (Aba) to (Alfm) show the development of a disturbance wave system thus 
generated. Only those waves generated when a blade of the rotor is aligned 
with a blade of the stator are portrayed >n these figures and then only a por- 
tion of each wave, extending a distance x , as shown in Figure (A4c) , are 

w 

shown for clarity. One complete revolution of the rotor, an angle of rotation 
of 2 tt, is represented over 12 intervals. The time increment for each interval 
? s 


2 it = 2tt 

= 12w 


(At) 


where is the number of blades in the stator, N 2 is the number of blades in 
the rotor and to is the rotational speed. The corresponding angular change for 
each At is 


A<}> 


2 * 

12 


(A2) 
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FIGURE a 3. ILLUSTRATION OF CYCLIC ALGORITHM FOR STAGE WITH UNEQUAL NUMBER OF 
BLADES IN STATOR AND ROTOR (N, = 3, N,*4) 



r 


Each time two blades are aligned, indicated by an arrow in Figure (A4a) 
through (A4m) , a pulse, shown as a large dot at x = 0 above the alignment, is 
created; this pulse then travels out to form a cylindrical disturbance wave 
with increasing radius in time. Figure (A4a) has a pulse labelled (1) cre- 
ated at t = 0 where the two reference blades (cross hatched) are aligned. 

This pulse forms the wave, labelled (l) shown in Figure (A4b) after 
At * 2Tr/(l2o>) or A<f> - 2ir/12, and a new pulse, labelled (2), is created where 
the next two blades are aligned (see arrow). These d ? sturbances , (T) and (2), 
travel out to the positions shown in Figure (A4c) and a new pulse, labelled 
(3), is created where the next set of blades are aligned, (see arrow) . The 
subsequent Figures (A4d) through (A^m) show the waves travelling out further 
and new ones being created similarly to Figures (A^a) through (A**c) . The 
first three waves are identified in Figures (AUa) through (AAd) and (A%i) . 
Figure (AAm) can be considered as an asymptotic picture where the wave front 
in the whole periphery at a distance of x = (a-u) • 1 2At is almost flat. If 
asymptotic configurations at three times are now considered, say at times 
t Q , t + At and t + 2At, as portrayed in Figures (A5a) through (A5c) , several 
important conclusions can be drawn. Consider first the geometry of the con- 
figurations. In Figure (A5a) the cross hatched reference or first set of 
blades are aligned (note the control volume indicated by the dot-dash border, 
to the top and right of it). In Figure (A5b) , representing a time At 
later, the next or second set of blades to the right are aligned (also note 
the corresponding control volume as described above). Finally, in Figure 
(A5c) 2A t after reference time t the third set of blades are aligned (this 
also corresponds to the two blades at the far left of Figure (A5c) , thus the 
control volume to the right of this set of blades shall be noted). Each of 
these dot-dash control volumes correspond to the same relative geometry and 
it can be seen by direct comparison that the asymptotic wave pattern is 
identical in each. The same thing is true at all other possible relative 
blade positions, (which implies all other times) that are encountered as 
the reference blade of the rotor moves through one complete blade spacing of 
the stator. It is thus possible to construct the entire asymptotic solution 
of the periphery at any one instant (and, therefore, at all instances) from 
all of the asymptotic periodic solutions at successive times found in one 
blade spacing control volume as the reference blade of the rotor travels 
through one blade spacing of the stator. 
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FIGURE AA. STARTING PROCESS SHOWING DEVELOPMENT OF DISTURBANCE 

WAVES FORMED WHEN BLADES ARE ALIGNED (ACTUAL SOLUTION) 
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Next consider the dash-dash, , control volume always lying to the 

right of the stator reference blade. This is the reference control volume 

within which all solutions for the flow upstream of the stator will be found 

as a function of time. Consider now the asymptotic wave pattern which would 

be found to the right of the reference control volume at time t^+2At in Figure 

(A5c) . It is obvious that it is identical to the wave pattern within the 

reference control volume in Figure (A5b) at time t Q +At or one At time step 

earlier (arrows “A"). Next consider the asymptotic wave pattern which would 

be found to the left of the reference control value at time t +2At i.e., in 

o 

Figure (A5c). Again it is obvious that it is identical to the wave pattern 
within the reference control volume in Figure (A5a) at time t or two At 
time steps earlier (arrows "B"). In view of this fact, the appropriate por- 
tions of the earlier solutions found in the reference control volume can be 
applied as boundary conditions to the reference control volume at the current 
time before proceeding with the computation for the next At. This is carried 
out, refer also to Figure (A3a) and (A3b) , by equating the values at the 
first row of exterior grid points, say line 5, to the values that existed At 
time step earlier at correspond ing points, in this case line (3» when the geo- 
metrical configurations with respect to lines a and B were the same. A similar 
procedure is carried out for lines a and y with its own phase lag. In addition, 
all other horizontal boundaries (those between the blade rows and downstream of 
the rotor) and vertical boundaries which require phase lags for proper speci- 
fication are specified analogously. These vertical boundaries are on the 
downstream side of domain k and the upstream side of domain 5. 

Consider now the starting process for the computational model. Figure 
(A6a) through (A6m) portrays the computational starting problem comparable to 
the actual physical one shown in Figure (A^a) through (A4m) . The successive 
application of the boundary conditions in this manner puts a numerical phase 
lag into the solution through the imposed boundary conditions as the computa- 
tion is started. As in the actual physical starting process discussed earlier, 
blade alignment , indicated by an arrow between blades, produces a pulse, in- 
dicated by a dot and identified in the first three waves in the Figures (A6a) 
through (A6d) and in Figure (A7m) . Initially there are no previously computed 
boundary conditions with phase lag to apply, hence, the "tagged 1 * disturbances 
created in the reference control volume are of minimal accuracy and as such 
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FIGURE A5 ASYMPTOTIC DISTURBANCE PATTERNS GENERATED 

WHEN A ROTOR BLADE IS ALIGNED WITH A STATOR 
BLADE. 
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FIGURE A6. STARTING PROCESS SHOWING WAVE DEVELOPMENT OF DISTURBANCE 
WAVES FORMED WHEN BLADES ARE ALIGNED (SOLUTION WITH 
BOUNDARY CONDITIONS HAVING PHASE LAG). 
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are indicated by using dotted curves. When a more accurate boundary condition 
is available, i'.e., when each side of the reference control volume can be 
specified from the solution at an earlier time, the accuracy of the solution 
within the control volumewill be improved. This Is indicated by dashed lines, 
as for example in Figures (A6e) and (A6f) etc. The boundary conditions with 
phase lag are applied as already discussed previously and the arrows in Figures 
(A6a) through (A6d) <^» ve. examp les of the phase lag for several cases. As the 
solution proceeds in time it also becomes more accurate due to the repeated 
application of increasingly accurate boundary conditions; this is indicated 
schematically by an increasing solidity of the lines depicting the wave fronts 
in Fi gure (A6) . 

Ultimately all solutions of lesser accuracy propagate out of the computa- 
tional domain and only the increasingly more accurate disturbances remain. The 
result is that the waves associated with the approximated starting process are 
lost and the solution becomes asymptotic in the same fashion as shown in Figures 
(A5 a) through (A5c) . The cri terion which determines the improvement is not 
the angular distance that the rotor has travelled, but is the number of times 
the boundary cond i ti ons ■ have been applied; consequently, the larger the number 
of blades, the better is the solution after a complete revolution of 2ir. An 
asymptotically periodic solution should be attained when successive solutions 
at the interior points possess periodic relationships and phase lags with re- 
spect to each other which precisely match those imposed by the application of 
boundary conditions at the virtual grid points. 

Horizontal Boundaries 

The "horizontal*' (i.e., streamwise) boundary points are specified by the 
technique discussed above which is based on the fact that these surfaces are 
"planes of equivalence" with a phase lag. Figure (A7) shows a typical domain. 

A row of "virtual" grid points in this case k = 2 and k = KS+2 is placed one 
mesh pbint outside the calculation domain; at these points the flow conditions 
are equated to those existing at corresponding grid points one mesh point in- 
side the opposite boundary, in this case K = KS-1 and K = 5 respectively, at a 
time when the blade positions at that opposite boundary were the same as cur- 
rently exist along the subject boundary. (The doub 1 e- va 1 ued boundary rows 
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FIGURE A7. SCHEMATIC OF A TYPICAL TRANSFORMED SQUARE DOMAIN f 
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K=3, 4 and K=KS, KS+1, have been introduced to account for jumps in proper- 
ties across the slipstreams.) 

Specification of the phase lags for these, hori zonta 1 boundaries is de- 
scribed as follows. Let there be N blades in the upstream blade row and M 
blades in the downstream blade row with N _< M <_ 3N/2. Also let n be the num- 
ber of time steps required for one blade spacing of the larger pitch to be 
traversed and m the number of time steps required for the smaller pitch to 
be traversed. Then: 

Nn = Mm 

There are four sets of boundaries to be considered: 

A. The upper boundary of domains 1, 2 and 4, adjacent to 
K=KS , KS+1. 

B. The upper boundary of domains 6 and 7> adjacent to 
K=KS , KS+1. 

C The lower boundary of domains 1, 2 and 4, adjacent to 
K=3, 4. 

D. The lower boundary of domains 6 and 7, adjacent to 
K=3, 4. 

The specific phase shifts are as follows: 

(1) Set A at K=KS+2 = Set C at K=5 (n-m) time steps ago 

(2) Set B at K=KS+2 = Set D at K=5 (n-m) time steps ago 

(3) Set C at K=2 = Set A at K=KS- 1 (m-(n-m)) time steps ago 

(4) Set D at K=2 = Set B at K=KS-1 (m) time steps ago 

The values computed at every time step at K=5 and KS-1 must, therefore, be 
stored for later use. Since a considerable number of time steps are involved, 
this cannot be done using core storage and is done instead by disk or tape 
storage and recall in a non-random method. 


Virtual grid rows K=1 and K=KS+3 are required in domains k, 6 and 7 to 
evaluate the spatial derivatives needed at K=2 and K=KS+2 in the characteristic 
solutions along the slipstreams. The procedure to define the solution at these 
exterior rows follows exactly the format outlined above, e.g.. Set A at 
K=(KS+3) is equated to Set C at K=6 (n-m) time steps ago. 

Entrance and Exit Boundaries 

Solutions at the virtual grid points on vertical boundaries J=1 and 
J=JS+1 in Figure (A7) (except, of course, the inlet and discharge stations 
which are true boundaries rather than virtual boundaries) are linearly inter- 
polated from the data in adjacent domains. Refer to Table Al for the following 
discussion for specific distribution of data at equivalent vertical boundaries. 
It is to be noted that in order to obtain spatially equivalent grid points for 
the vertical boundaries between domains k and 5 account must be taken of the 
relative motion of the domains; interpolation of data as well as application 
of phase shifts are required. 

In the case of equal numbers of blades in adjacent rows, matching flow 
variables along the interface between domains k and 5 is straightforward; in 
any other case certain subtleties occur in the process of minimizing the com- 
puter storage and time requirements. Again a column of grid points external to 
domains k and 5 are considered. In the regions of overlap shown in Figures 
(A8a) and (A8b) , it is clear that the exterior column of points from domain 
overlaps the interior of domain 5, and vice versa. However, for the regions 
of "non-overlap" determination of the flow variables at these exterior grid 
points is accomplished by a phase lag technique, analogous to the periodic 
boundary technique discussed above, which is based on the blade positions during 
a full cycle of movement of domain 5 relative to domain b. In the case of 
equal spacing the time-delay is zero, and the regions of non-overlap labelled 
"left" and "right" in Figure (A8a) have a direct correspondence. On the other 
hand, in the case of fractional (non- integer) spacing the matching of conditions 
along the region of non-overlap must be based on a correspondence of relative 
blade positions, which introduces a phase lag in the application of data on the 
interface between domains ^ and 5* 
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This phase is analogous to that discussed earlier for the horizontal 
boundaries. The only added feature is that up to three separate sets of 
data each with its appropriate phase shift must be applied since the overlap 
can be such as shown in Figure (A8c) which clearly shows three regions where 
data must be provided. Thus the J=JS values for domain 4, refer to Figure 
(A7) , are equal to the J=2 values of domain 5 at times having up to three 
different phase shifts, as discussed earlier, and the J=1 values for domain 
5 are equal to the J=JS-1 values for domain 4 at times having up to two dif- 
ferent phase shifts. Even though only two are needed, three are specified 
since the non-overlap region can be either above or below the larger blade 
gap. The interpolation merely ignores the extra data. 
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TABLE Al 


EQUIVALENT VERTICAL BOUNDARIES, WITH AND WITHOUT PHASE LAG 


DOMAIN 


DOMAIN 

AND J INDEX 


AND J INDEX 

(Exterior) 


( Interior) 

1 (JS) 

= ' 

2(2) 

2(1) 

= , 

Interpolated between 
1 (JS- 1 ) and 1 (JS) 

2 (JS) 

= 

3(2) 

3 (D 

= 

2 (JS- 1 ) 

3 (JS+1 ) 

= 

f Interpolated from Domain 4 

l 6(3) 

4 (2) 

= 

3 (JS) 

4 (JS) 

= 

5(2) 

5 (1) 


Interpolated from Domain 4 

5 (JS+1) 

= 

6(3) 

6 (2) 

= 

5 (JS) 

6 (JS+1) 

= 

Interpolated between 
7(2) and 7(3) 

7 (2) 

= 

6 (JS) 


COMMENTS 


Only if Domain 1 is active 
Only if Domain 1 is active 


Two blade rows 
One blade row 

Only if Domain 4 is active 
(i.e., two blade rows). 


Interpolated with or without 
phase lag (IBLEQ^O or l), 
only if Domain 4 is active. 


Interpolated with or without 
phase lag (lBLEQ=0 or l), 
only if Domain 5 is active. 

Only if Domain 5 is active 
(i.e., two blade rows) 


Only if Domain 7 is active 


Only if Domain 7 is active 
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